亞臨界區圓柱繞流相干結構壁面?；旌蟁ANS/LES 模型

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

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

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

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 理論模型

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

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

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

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

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

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

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

對DES 模型[29],定義為

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

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

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

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

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

其中,fd為一個經驗性混合函數,具體形式如下:

其中,dw為網格計算點到壁面的距離,κ為馮卡門常數,取值為0.41.

函數fd的分布特點為,在近壁層的某個范圍內(與網格及當地流場有關)等于0,即在此區域內混合長度=lRANS,此時DDES 模型還原為RANS模式.在此區域外的計算區域,混合長度決定于兩個尺度lRANS和lLES的相對大小,此時DDES 與DES的表現是一樣的.

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

hmax為計算單元的最大網格步長.

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

此外,由文獻[35]可得

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

由(16)式可得

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

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

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

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

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

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

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

其中,Cw=0.15,hwn為壁面法向網格步長.

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

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

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

克服IDDES 模型缺陷的一種有效途徑是使其具有壁面?；鬁u模擬的能力.有鑒于此,本文構造一種新的混合函數fnd,使新所構造的混合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 不穩定性結構這類特殊流動現象,其對湍動能的解析能力至少需要達到75%[2],即rk≤0.25 .

定義kc為截斷波數,它表示LES 的濾波寬度,可由當地網格尺寸Δ*確定如下:

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

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

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

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

對當地網格尺度Δ*,一種合理的選擇[39]為Δ*=hmax.此時,rk可以改寫為

根據Nyquist-Shannon 樣本定理,為了能夠準確解析相關尺度的湍流結構,湍流含能尺度Ltur需要滿足如下條件[39]:

其中,Nr為某個長度尺度結構能夠被解析的網格單元數量,它與單位波長的波能夠被重構所需的樣本點數量一致.

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

有鑒于此,設rk1≤.本文將當地網格滿足條件rk＞rk1的計算區域定義為WM-HRL 模型的完全RANS 區域,即

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

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

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

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

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

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

為此,下面繼續構造rk的一種僅依賴于當地網格尺寸的計算公式.Pope[40]指出,在對數律層區,Ltur與dw成正比關系,即Ltur=Cwdw.在高雷諾數的情況,Han等[41]指出,Cw可近似取為Cw≈2.5 .由此可見,(30)式可以改寫為:

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

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

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

其次,當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 混合轉換的區域定義為

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

根據如上所述兩個nrk1-Q和nrk2-Q準則,可構造一種新的混合函數fnd如下:

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

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

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

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

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

下面證明,在(42)式的條件下,由(41)式構造的函數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 模型為所謂的代數壁面?；鬁u模擬(Alg WMLES)[34],其基本思想是對Smagorinsky 型LES 的亞格子渦黏系數乘以一個阻尼系數Fr,即:

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

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

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

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

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

3 計算模型與設置

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

圖1 計算區域設置Fig.1.Computational domain schematic.

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

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

圖2 計算網格剖面Fig.2.Computational grid configuration.

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

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

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

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

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

表2 文獻中雷諾數Re=3900 下圓柱繞流場相關統計量的實驗和數值結果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 模型,濾波控制函數fk定義為fk=ku/k.Lehmkuhl等[10]采用了兩種DNS模式.其中,“Mode H”表示“高能模態”,利用該模態的DNS 得到V 形結果.而“Mode L”表示“低能模態”,利用該模態的DNS 得到U 形結果.

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

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

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

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

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

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

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

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

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

表3 監測點坐標信息Table 3.Coordinate information of the probes.

圖3 剪切層小尺度K-H 不穩定性結構監測點分布Fig.3.Location configuration of the probes for the small scale K-H instability structure in the shear layer.

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

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

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

4 數值結果與分析

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

4.1 WM-HRL 模型參數影響規律

首先討論LES 啟動邊界ΓLES位于剪切層小尺度K-H 不穩定性結構發生區域的情況,并考慮兩種情況:第一種情況為=105.8,rk2=0.1556 ;第 二種情況為=113.9,rk2=0.1484 .在每種情況下,均設置11 種RANS 結束邊界ΓRANS的位置及其相應的rk1,如表4 所示.其中,=7.9位于黏性底層,=13.3,14.9,18.4,20.4,27.1 位于過渡層,=38.4,41.7,49.2,72.7 位于對數律層但在剪切層K-H 不穩定區外,而=91.2則位于剪切層K-H 不穩定性結構發生區域內.

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

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

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

表6 當 ΓLES 位于過渡層時,相關流場統計量的數值結果Table 6.Numerical results for flow statistic characteristics when ΓLES is located in the buffer layer.

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

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

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

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

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

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

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

4.2 一階和二階統計量特性

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

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

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

4.2.1 一階統計量特性

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

圖4 圓柱表面周向壓力系數 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 模型所得沿圓柱表面周向壓力系數Cp分布的數值結果之間的差異均很小.對圖4(a)的情況,本文計算結果與Norberg[43]實驗結果相比的相對誤差最大為10.3%,而Lehmkuhl等[10]利用DNS(Mode L)計算結果與Norberg[43]實驗結果相比的相對誤差最大為3.1%.對圖4(b)的情況,本文計算結果與Norberg[43]實驗結果相比的相對誤差最大為14.0%,而Lehmkuhl等[10]利用DNS(Mode H)計算結果Norberg[43]實驗結果相比的相對誤差最大為0.6%.

在圖5 中,分別給出了前述8 種工況下沿尾流中心線平均流向速度剖面特性的數值結果.由圖5可知,在圓柱表面正后方中心線上的點P(0.5D,0)處,平均流向速度的值為0,在達到一個負的最小值后開始增大,并在中心線上的點Q(Lr+0.5D,0)處又變為0.其中,Lr/D為回流區長度.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)所示.此即在雷諾數Re=3900 下圓柱繞流兩個實驗及DNS 計算中出現兩類不同回流區分支結構的情況.

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

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

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

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

圖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 模型所得該站位處平均流向速度分布的數值結果之間的差異均很小.在位于剪切層K-H 不穩定性結構發生區域的站位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 位于回流區中.由圖6可見,Parnaudeau等[18]和Lourenco 和Shih[27]的實驗結果均獲得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)處平均橫向速度剖面特性的數值結果.

圖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 模型所得圓柱后方不同站位處平均橫向速度分布的數值結果差異均很小.

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

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

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

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

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

4.2.2 二階統計量特性

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

圖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 不穩定性結構發生區域的站位x1/D=1.06 處,本文計算結果與Parnaudeau等[18]的實驗結果及Lehmkuhl等[10]利用Mode L 之DNS 計算結果相比,主要差異在兩個“峰”處.在位于回流區的站位x1/D=1.54,2.02處,本文計算結果與Parnaudeau等[18]的實驗結果及Lehmkuhl等[10]利用Mode L 的DNS 計算結果相比,在兩個“峰”處的值小于相應的實驗值及DNS 計算值.

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

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

圖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 不穩定區的站位x1/D=1.06 處,本文計算結果與Parnaudeau等[18]實驗結果相比的最大相對誤差分別為23.6%,7.7%,15.7%,9.1%,而Lehmkuhl等[10]利用Mode L 之DNS 計算結果與Parnaudeau等[18]實驗結果相比的最大相對誤差為14.5%.

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

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

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

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

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

圖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 不穩定性發生區域的站位x1/D=1.06 處,本文在Case AU—Case DU 這4 種工況下的計算結果均與Parnaudeau等[18]的實驗結果具有良好的一致性.在位于回流區的站位x1/D=1.54 處,本文對Case AU 和Case BU 的計算結果與Parnaudeau等[18]的實驗結果也具有良好的一致性,但在 |x2/D|＜0.48 的范圍內本文對Case CU和Case DU 的計算結果,與Parnaudeau等[18]的實驗結果有一定的差異,最大相對誤差分別為38.7%和14.8%.在位于回流區的站位x1/D=2.02處,本文對Case BU—Case DU 的計算 結果 與Parnaudeau等[18]的實驗結果具有良好的一致性,但在 |x2/D|＜0.48 的范圍內本文對Case AU 的計算結果,與Parnaudeau等[18]的實驗結果有一定的差異,最大相對誤差為30.8%.

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

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

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

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

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

研究發現,對本文所構建的WM-HRL 模型,在RANS/LES 混合轉換區邊界位置及其兩個指標參數nrk1-Q和nrk2-Q準則取值的某些組合下,計算所得二階統計量與相關文獻的實驗結果相比,存在較大的誤差,其主要原因如下.

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

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

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

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

4.3 相干結構頻譜及流場特征

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

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

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

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

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

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

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

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

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

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

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

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

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

最后討論圓柱繞流場中兩類相干結構的流場特征,計算工況選取Case AU—DU 及Case AV—DV 兩類情況,結果如圖14 和圖15 所示.其中,圖14 的左列為x3/D=π/2 斷面上展向渦量ω3=?u2/?x1-?u1/?x2的云圖,右列為流向速度云圖,且圖中白色垂直虛線為回流區末端位置(即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 (左)均清晰地可見附著于圓柱上下側壁的兩個層流分離剪切層結構,表現為兩條位置狹小的窄帶結構.同時,從圖14 (右)則可清晰地觀察到圓柱尾部因流動分離而形成的兩個回流區分支結構.具體來說,第1 個分支結構發生在x1/D=1.06 處平均流向速度剖面為U 形的情況,回流區長度較長;而第2 個分支結構發生在x1/D=1.06 處平均流向速度剖面為V 形的情況,回流區長度較短.

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

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

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

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

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

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

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

5 結論

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

通過構造一個僅與當地網格空間分布尺寸相關的湍動能解析度指標參數rk,提出了確定該WM-HRL 模型兩個邊界位置和3 個區域湍動能解析度的nrk1-Q和nrk2-Q準則.其中,nrk1-Q準則用于確定WM-HRL 模型中RANS 模式結束邊界位置,及其在RANS 區湍動能的最大解析度(即rk1),而nrk2-Q準則用于確定WM-HRL 模型中LES 模式啟動邊界位置及其在LES 區對湍動能的最小解析度(即rk2).

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

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

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

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

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

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

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

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

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

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