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

1) (上海交通大學,海洋工程國家重點實驗室,上海 200240)

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

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 理論模型

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

其中,頂標“—”表示雷諾時均.τki為Reynolds 應力,采用如下Boussinesq 近似進行計算:

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

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

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

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

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

對DES 模型[29],定義為

其中,lRANS為RANS尺度,lLES為LES尺度,可分別表示為

其中,CDES為模型常數(shù),Δ為LES 濾波寬度.

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

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

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

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

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

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

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

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

對DES,DDES 和IDDES 模型,當其進入LES后,在局部平衡流條件下,SST 模型k方程中的生成項與色散項相等[35],即

此外,由文獻[35]可得

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

由(16)式可得

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

本文對雷諾數(shù)Re=3900 下圓柱繞流場問題,采用WM-HRL 進行數(shù)值模擬與分析.其中,Re=U∞D(zhuǎn)/v,U∞為來流速度,D為圓柱直徑.計算區(qū)域為矩形區(qū)域,如圖1 所示.具體設(shè)置如下:圓柱底面中心位于坐標原點 (0,0,0),入口邊界位于x1/D=-10 ;出口邊界位于x1/D=15 ;前后兩個邊界分別位于x2/D=±4 ;展向高度為L3/D=π .

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

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

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

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

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

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

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

表1 在雷諾數(shù)Re=3900 下圓柱繞流文獻中所用計算模型與網(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)格分辨率略低于前三篇文獻的情況,但水平方向的網(wǎng)格分辨率高于Breuer[15].Pereira等[2]和D'Alessandro等[30]采用的展向網(wǎng)格分辨率最低,兩者的總網(wǎng)格量接近.本文所采用的展向網(wǎng)格分辨率與Tremblay[8]和Breuer[15]的一致,但水平方向的網(wǎng)格分辨率均要低于這兩篇文獻的情況.總之,本文所用計算網(wǎng)格數(shù)量均少于DNS[10],LES[8,15],PANS[2,24]及DES類模型[24,30]的計算網(wǎng)格數(shù)量.

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

表2 文獻中雷諾數(shù)Re=3900 下圓柱繞流場相關(guān)統(tǒng)計量的實驗和數(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 中,對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 可知,對Pereira等[2](PANS) (fk=0.25,0.5)的情況,與表中所列實驗結(jié)果[18,27]相比,分離角的計算結(jié)果均明顯偏小.對Luo 等(PANS)[24](fk=0.5)的情況,分離角?s、阻力系數(shù)Cd及基線壓力系數(shù) -Cpb的計算結(jié)果與表中所列實驗結(jié)果結(jié)果[27]相比均明顯偏大,相對誤差分別達到9.176%,37.76%和63.333%.對D'Alessandro 等(SA-DES)[30]的情況,阻力系數(shù)Cd的計算結(jié)果與表中所列實驗結(jié)果[27]相比均明顯偏大,相對誤差達22.7%.

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

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

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

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

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

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

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

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

表3 監(jiān)測點坐標信息Table 3.Coordinate information of the probes.

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

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

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

在圖2 所述網(wǎng)格系統(tǒng)設(shè)置下,利用(36)式進行計算,結(jié)果表明:在對數(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 模型對該剪切層小尺度K-H 不穩(wěn)定結(jié)構(gòu)為完全LES 解析模式,且其對湍動能的解析度至少為80%.

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

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

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

首先討論LES 啟動邊界Γ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的位置及其相應的rk1,如表4 所示.其中,=7.9位于黏性底層,=13.3,14.9,18.4,20.4,27.1 位于過渡層,=38.4,41.7,49.2,72.7 位于對數(shù)律層但在剪切層K-H 不穩(wěn)定區(qū)外,而=91.2則位于剪切層K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域內(nèi).

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

表5 當 ΓLES 位于剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域外且在對數(shù)律層內(nèi)時,相關(guān)流場統(tǒng)計量的數(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 可知,只有當RANS 結(jié)束邊界ΓRANS位于過渡層,且rk1≤=0.63 時,才能同時準確獲得和Lr/D的計算結(jié)果.在此條件下,當計算得到V 形速度剖面時,Lr/D的計算結(jié)果與Lourenco和Shih[27]的實驗結(jié)果相比,相對誤差最大為-9.32%,而的計算結(jié)果與Lehmkuhl等[10]的DNS 之Mode H 計算結(jié)果相比,相對誤 差最大為12.68%.當計算得到U 形速度剖面時,Lr/D的計算結(jié)果與Parnaudeau等[18]的實驗結(jié)果相比,相對誤差最大為-12.58%,而的計算結(jié)果與Lehmkuhl等[10]的DNS 之Mode H 計算結(jié)果相比,相對誤差最大為-20.15%.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖7 分別給出了前述8 種工況下在3 個不同站位(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 可見,無論是對Case AU—Case DU 這4 種工況,還是對Case AV—Case DV 這4 種工況,利用本文WM-HRL 模型所得圓柱后方不同站位處平均橫向速度分布的數(shù)值結(jié)果差異均很小.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖11 在Case CU 工況下,在P1—P12 監(jiān)測點處橫向脈動速度的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)測點處橫向脈動速度的Lomb譜Fig.12.Lomb spectrums of the cross-stream fluctuation velocities at different probes P1-P12 for the Case CV.

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

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

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

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

圖13 在Case CU (第1 和第2 列)和Case CV (第3 和第4 列)工況下,在P13—P18 監(jiān)測點處流向(第1 和第3 列)及橫向(第2 和第4 列)脈動速度的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 可見,在Case CU 和Case CV 兩種工況下,對P13—P18 這6 個監(jiān)測點,無論是在其流向脈動速度的Lomb 譜中,還是在其橫向脈動速度的Lomb 譜中,均清晰地出現(xiàn)兩個頻率較低的譜峰,一個譜峰對應渦泄頻率(),另一個譜峰對應渦泄頻率的倍頻().

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

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

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

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

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

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

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

圖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 (左)均清晰地可見附著于圓柱上下側(cè)壁的兩個層流分離剪切層結(jié)構(gòu),表現(xiàn)為兩條位置狹小的窄帶結(jié)構(gòu).同時,從圖14 (右)則可清晰地觀察到圓柱尾部因流動分離而形成的兩個回流區(qū)分支結(jié)構(gòu).具體來說,第1 個分支結(jié)構(gòu)發(fā)生在x1/D=1.06 處平均流向速度剖面為U 形的情況,回流區(qū)長度較長;而第2 個分支結(jié)構(gòu)發(fā)生在x1/D=1.06 處平均流向速度剖面為V 形的情況,回流區(qū)長度較短.

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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

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

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

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

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

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

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

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