基于噴流擬序結構預測的SGS 模型比較研究1)

劉琪麟 賴煥新

(華東理工大學機械與動力工程學院,上海 200237)

引言

擬序結構普遍存在于湍流中[1-3],并影響湍流的發生、輸運與耗散等過程,對湍流擬序結構的研究、對于理解湍流機理、分析流動相互作用有重要的意義.從大渦模擬來看,正確預測湍流擬序結構的關鍵就在于算法和亞網格尺度(sub-grid-scale,SGS)模型的準確性.

在算法方面,傳統的大渦模擬方法求解網格濾波的控制方程,并采用SGS 模型模化亞網格尺度流動與網格尺度流動之間的相互作用.雖然近年來部分學者利用數值離散的特點,也發展出不需要添加模型的隱式大渦模擬方法[4-5],但為了針對性地討論應用中亞網格模型的選擇問題,本文仍采用高精度有限差分結合亞網格模型的大渦模擬經典算法,并討論SGS 模型的性能.對大渦模擬算法更全面的介紹可以參考Sagaut 等[6]的綜述.

目前文獻中廣泛出現的亞網格尺度模型是經典的常系數Smagorinsky 模型(Smagorinsky model,SM),已有眾多文獻報道該模型過度耗散了脈動動能[7-9].另一類是基于Germano 恒等式和Lilly 最小二乘法來動態調整模型系數的動態Smagorinsky 模型(dynamic Smagorinsky model,DSM),這類模型的局限性在于理論上要求在流場均勻的方向上計算平均的模型系數,而這一平均過程缺少物理上合理的解釋,且造成分布式并行計算中不同處理器之間數據的頻繁調用與交換,降低計算效率.從計算的角度,較為理想的模型是局部模型,計算效率高且能較好地克服SM 對湍流黏性估計過高的缺點.Piomelli和Liu[9]提出了局部動態模型(localized dynamic Smagorinsky model,LDSM),去除了DSM 對流動中需存在各向同性方向的要求,而是利用流動的歷史信息來避免平均過程,計算效率得到提高.Lenormand等[10]提出的選擇多尺度模型(selective mixed-scale model,SMSM),在邊界層流動[10]、翼型繞流[11]和空腔流動[12]的預測中具有良好應用.這一模型借助渦矢量方向的變化來判斷流動的狀態是否發生了轉捩,同時反映湍流的間歇性.Kobayashi 等[13-14]則利用流場中黏性耗散與擬序結構的關聯性,提出了擬序結構模型(coherent-structure Smagorinsky model,CSM),并應用于各向同性湍流、湍流通道流以及橫向射流的計算中,其模型簡單且具有類似于動態模型的性能.Koyabashi[13]同時改進動能模型,提出擬序結構動能模型(coherent-structure kinetic-energy model,CKM).這些局部模型都有望克服傳統模型的不足,但對于預測可壓縮噴流中的擬序結構,目前尚未看到有關它們預測性能的分析與評估的文獻.

典型的噴流擬序結構有流向渦、高低速流體間的條帶[15]、剪切層邊沿的卷吸流動結構[16]等,表現為長時間保持一定形狀的大尺度質量團,且在時間和空間上具有強相關性[17].Lumley[18]最早提出采用本征正交分解(proper orthogonal decomposition,POD)方法提取擬序結構,且在自由剪切流、分離流、邊界層轉捩以及湍流通道流中得到廣泛的應用.POD 方法的目標是針對物理場尋找一組空間正交的基函數,并按照基函數對物理量方差的貢獻排序,例如對于脈動速度場的研究就是基函數對湍動能的貢獻[19-20].本征正交分解的特點使得只需要少量的POD 模態就能夠提取出流場中能量占主導地位的流動模式.此外,POD 分解的基函數直接由湍流場確定,不同于傅里葉變換或者小波變換需要人為指定基函數,因此表征的擬序結構更為客觀[17].在POD 算法方面,雖然快照POD 算法能夠顯著減少計算量,但是有部分高階模態沒有參與計算.基于矩陣奇異值分解(singular value decomposition,SVD)的經典POD 算法包含完整的模態信息,比快照POD算法對舍入誤差的魯棒性高[19-20].此外數學軟件Matlab 中提供了SVD 的庫函數,也便于編制后處理程序.

本文對亞聲速可壓縮平行噴流進行大渦模擬.在驗證計算結果的基礎上,研究SM,CKM,SMSM,LDSM 和CSM 模型的特性,分析各個模型預測的瞬時渦結構.同時采用基于SVD 的POD 分解方法提取噴流擬序結構并分析各模型預測結果的差異,為可壓縮噴流的大渦模擬提供參考.

1 數值計算與方法

1.1 亞網格尺度模型

控制方程組為Favre 質量加權濾波的三維非穩態可壓縮N-S 方程組[21].方程中的亞網格尺度未知量由亞網格模型給出,包括亞網格尺度應力τij

以及亞網格尺度熱通量

其中ρ,T,νsgs以及ui(i=1,2,3) 分別為密度、溫度、模化的亞網格黏性系數以及速度矢量的3 個笛卡爾分量.γ 是量熱完全氣體的比熱比,空氣取γ=1.4.M為馬赫數,湍流普朗特數常數Prt=1[22],而,頂標“-”表示網格濾波,“～”表示Favre 質量加權的網格濾波,所有計算變量均為網格濾波后的變量,因此除非特殊說明以下省去頂部符號.

Smagorinsky 模型(SM)中的亞網格尺度應力為

δij=0.對于噴流,取模型系數C1=0.17[23].τkk的模化根據文獻[7]的建議得出.其中=Li j?δi jLkk/3,頂標“^”表示測試濾波運算,Ct為上一時刻的模型系數,初始值可以根據SM 模型取Ct=C1,當前時刻的模型系數C2是未知量,且式(5)對C2是超定的.C2的估計值由式(5)兩邊取αij的縮并得出[9],如式(6)所示.LDSM 模型的亞網格尺度應力與亞網格黏性系數通過將式(4)中的替換為得出

選擇多尺度模型(slective mixed-scale model,SMSM)是根據局部渦矢量方向的脈動來調整亞網格黏性系數

其中C3=0.06,α 取0.5[24].為了反映湍流的多尺度特性,引入小尺度脈動的動能=(1/2),其中脈動速度=?.進一步,對網格濾波尺度的瞬時渦量進行測試濾波,得出.和兩個向量之間的夾角 θ 反映了渦運動的脈動強弱.對于均勻各向同性湍流,在夾角 θ 的概率密度函數的峰值位置處有θ0=20°[24],據此,可判斷局部流動的狀態.模型引入了選擇函數fθ0,當 θ> θ0時,取fθ0=1,表示渦方向在兩個濾波尺度上變化較大,是無規則的完全湍流狀態;否則,當 θ<θ0時,取fθ0=tan(θ/2)/tan(θ0/2),意味著湍流逐步減弱,亞網格黏性也應相應減弱,直到在層流區僅保留分子黏性.因此,SMSM 模型的亞網格黏性系數為

擬序結構模型是利用擬序結構與湍流中耗散的相關性來調整亞網格黏性系數vsgs=CΔ2|S(u)|,模型系數C=C4|FCS|3/2中常數C4=0.05[13],FCS=Q/E為擬序結構函數，其中

為速度導數Dij=?ui/?xj的第二不變量,它主要由Di j中反映當地旋轉的反對稱部分Aij和反映當地變形的對稱部分Bi j組成,而為Di j的范數.其中

基于直接數值模擬(direct numerical simulation,DNS),流場的相關性分析表明,Q的絕對值與擬序小尺度渦附近能量耗散的強弱存在關聯[13],因此|FCS|反映了流場中黏性與耗散的局部變化,而|FCS|的冪次“3/2”是根據不可壓縮壁面流動中C∝y3,Q∝y2以及E∝ 常數(y是壁面的法向)而確定的.因此,CSM 模型的亞網格黏性系數為

擬序結構動能模型是在擬序結構模型的基礎上,在亞網格黏性系數中顯式地包含亞網格尺度動能ksgs,并做量綱調整后得到

至此,可以通過式(2)、式(4)來計算亞網格尺度熱通量和亞網格尺度應力.雖然文獻[12]指出亞網格尺度應力中各向同性的部分 τkk相對于熱力學壓力可以忽略,但根據Vreman[7]的建議,在基于SM 的LDSM 以及CSM 模型中仍然模化 τkk.

1.2 計算對象與數值離散

本文采用Hu 等[22]的亞聲速平行噴流為考核算例.噴流出口雷諾數ReJ=DJρ∞UJ/μ∞=2000,馬赫數MJ=UJ/c∞=0.9,其中DJ=1、UJ=1 分別為歸一化的噴口寬度和噴口速度,參考值ρ∞,μ∞和c∞分別為來流密度、動力黏度以及環境聲速.計算域與網格劃分如圖1 所示,流向(x)、橫向(y)與展向(z)的尺寸分別為Lx×Ly×Lz=24×15×3,采用笛卡爾網格進行離散,對應方向的結點數分別為 181×181×16.為了提高近噴口區域和噴流唇線附近的網格分辨率,采用雙曲正切函數對x和y方向的網格進行拉伸,而噴流在z方向采用周期邊界條件,劃分均勻的網格.基于Celik 等[25]提出的大渦模擬分辨率指標,由于課題組前期工作已證明該網格能夠分辨湍動噴流中75%以上的湍動能,考慮到噴流的雷諾數較低,該計算網格可以滿足大渦模擬的網格質量要求[26].

圖1 計算域與網格(3 個方向都每隔3 條網格線畫一條線)Fig.1 Computional domain and the Grid(One in every three lines is shown in each direction)

采用有限差分方法開展數值模擬,本文的計算源程序以Sandham 課題組[27-28]的(shock-boundary layer interaction,SBLI)代碼為基礎,計算精度和可靠性經過了充分的考核驗證[22,29].在此基礎上,本文進一步采用OpenMP 并行化,使該代碼可用于內存共享式的小型計算工作站.

空間離散的1 階與2 階導數均采用5 點4 階中心差分格式計算,邊界點采用Carpenter 提出的單側4 階穩定差分格式[30],從而保持空間離散的整體高精度.為了抑制差分計算引起的非物理振蕩,算法上采用了熵分裂方法,通過將對流項分裂為守恒形式與非守恒形式的加權求和來實現[31].同時,在x和y方向的計算邊界上,采用無反射邊界條件來抑制非物理振蕩波在計算域邊界上的反射,算法上通過保留走出計算域的特征波而將走入計算域的特征波幅值置零來實現[32].在z方向上采用周期邊界條件來模擬平行噴流在展向的無限延伸.

時間離散采用低存儲顯式3 階龍格庫塔算法,這一算法完成3 個子時間步的推進僅需要兩個存儲流場的數組,對內存資源的消耗低[21].為便于模型特性的研究比較,時間步長取定值Δt=0.003 075,且滿足數值穩定性條件CFL<1.在初始條件方面,給定噴口平面上流向速度u、溫度T以及壓強p的形線,并用噴口平面上的物理量來初始化全場.噴口的速度分布采用雙曲正切函數給定,并施加伴流速度U0=0.1,可得

在數值結果方面,雖然DNS 能夠準確給出流場結構的信息,但是DNS 對計算資源的消耗量高,而且流場后處理的數據量大.以文獻[22-33]中提供的DNS 數據為基礎,課題組前期工作已開展了計算代碼的考核[34-35]、LES 的網格分辨率分析[26]以及局部模型的研究[21,36],算法的可靠性與模型都已得到驗證.本文在已有工作的基礎上,進一步分析局部模型揭示的噴流擬序結構.

1.3 POD 方法

POD 方法采用一組正交基函數 ψ(x) 的線性組合來表示物理場 φ(x,t),組合系數為c(t),即

所得的基函數也稱為POD 模態.POD 分解的目的是使這種表示的誤差在平方平均的意義下最小,并且在給定誤差下最小化表示所需要的模態個數,從而得到最佳的基函數[20].這組最佳基函數能夠最優地表征原物理量的方差,對于脈動速度場，這種最優是模態反映的湍動能隨模態階數增大而迅速的減小[19],脈動速度的第一階POD 模態對湍動能的貢獻最大.詳細的理論介紹可參考文獻[20].

在計算中,可對流場的速度分量進行時間采樣,得到關于標量 φ 的N個流場樣本,對于網格點數為M=Mx×My的二維樣本,可將每個樣本都寫成M×1的列向量,由此得到M×N的樣本矩陣

在離散求解時,尋找最佳基函數的過程在數學上等價于求樣本矩陣F的SVD 問題[20].采用Matlab 中低存儲的“svd()”庫函數計算F的奇異值分解

其中,左矩陣 Ψ 的各列是正交的,Ψ=(ψ1,ψ2,···,ψN),它的第i列 ψi即為物理場的第i階POD 模態,右矩陣C的各列也是正交的,C=(c1,c2,···,cN),它的第i列ci即為模態 ψi對應的系數,反映POD 模態與時間有關的特征[19],而S為奇異值矩陣,S=dig(s1,s2,···,sN),S的對角元為從大到小排列的奇異值si(i=1,2,···,N),其他矩陣元為零,λi=代表了與第i階POD 模態對應的能量,且對各階模態可以定義能量貢獻率

以及累計能量貢獻率

另外,由于各階POD 模態 ψi是相互正交的,所以λ衰減越快,當前模態所含有的能量相比下一階模態的差就越大,當前模態所表征的物理場的相關性就越強.

2 SGS 模型的特性

2.1 平均流場

在分析湍流統計特征時,為排除初場對統計結果的影響,取無量綱時間t=48～ 216 的流場數據,對應流體以噴口速度流過2～ 9 倍的計算域長度,以計算時間步長Δt為采樣間隔,得到54 640 個流場樣本.圖2 以LDSM 模型為例,顯示了在唇線上無量綱位置x=2,6,8,10,20 處的流向脈動速度的功率譜,譜線具有一段固定的斜率且能觀察到慣性區,反映了湍流的能量級串現象,也說明流動已經達到了完全湍動狀態.

圖2 LDSM 模型預測的脈動速度無量綱功率譜Fig.2 Dimensionless power spectra of velocity signals predicted by the LDSM

噴流主流流體與環境流體存在速度差,形成剪切層.圖3 所示為噴流剪切層的發展.考慮噴流中心線速度Uc與伴流速度U0,圖3(a)采用速度半值寬δ0.5,即流向速度u= 0.5(Uc+U0) 時的橫向位置y,顯示出剪切層的發展呈現兩個線性階段,即初期的緩慢增長段與下游的快速增長段,這對應噴流勢流核末端x≈6 處的流動轉捩現象.對比5 種模型的結果,SM 模型預測的剪切層橫向發展程度較其他模型大,表明該模型過度耗散的特點.圖3(b)以LDSM 模型為例,在相似坐標中繪出了噴流不同流向位置的速度形線,相似坐標的橫向位置采用速度半值寬歸一化,平均噴流速度差〈u〉?Uc采用主、伴流的速度差ΔUc=Uc?U0歸一化,〈·〉 表示時間平均.圖3(b)中,速度形線在入口呈平頂狀分布,而在下游則快速發展為自相似的速度形線,可以看到x=4 時,自相似形線還沒有建立,但當x≥8 時,速度形線與實驗噴流[37-39]的自相似數據點吻合良好,也驗證了計算結果.

圖3 噴流剪切層的發展Fig.3 Developments of the jet shear layer

噴流流體與環境流體的強烈摻混伴隨著能量的耗散.對比5 種SGS 模型的亞網格尺度平均黏性耗散率,定義為.圖4 為在展向中心平面(z=1.5)上的分布.由于噴流入口施加的平均來流條件,流體噴出后在近噴口區基本為層流狀態,隨后在勢流核末端(x≈6)附近發生流動轉捩,并在下游發展為充分湍流.圖2 中,近噴口的脈動水平比勢流核及其下游的區域低1～ 2 個數量級,進一步表明來流的層流特征,此處與湍流脈動相關的亞網格尺度黏性耗散應趨近于零.但在圖4(a)中,SM 模型預測的在近噴口區內呈兩條平行的條帶,表明SM 模型對層流擾動敏感,不能正確反映湍流黏性耗散率的分布區域.圖4(b)～ 圖4(e)分別為CKM,SMSM,LDSM 以及CSM 的預測結果,這4 種局部模型雖然能夠正確地反映近噴口區的分布特征,在轉捩區(x=6～9)內呈現一對峰值,但是以CSM 模型為參考,CKM,LDSM,SMSM 以及SM 模型預測的峰值耗散不同,分別為CSM 模型預測值的3.0,1.8,1.3 和0.9 倍.其中CKM 預測的峰值接近其他模型的1～ 2 倍,稍后分析的分布對轉捩區內渦結構的影響.

圖4 對比在中心平面(z=1.5)上的分布Fig.4 Comparisons of the distributi on of ?in the center plane of the jet (z=1.5)

2.2 瞬時流場的渦結構

分析無量綱時間t=207 的瞬時流場,此時噴流以噴口速度流過8 倍以上的計算域長度,已經充分排除初場的影響.參考Hu 等[22]的部分DNS 結果,以噴流中心線速度〈u〉 衰減為0.99 時的流向位置得出平均的勢流核長度為x=6.43.圖5 為Q=0.1 的等值面,圖5 中紅色切片位于Hu 等DNS 結果[22]的勢流核末端.在預測流動轉捩方面,圖5(c)～ 圖5(e)中,SMSM,LDSM 以及CSM 模型都能預測大尺度展向渦在平均勢流核末端附近的破碎,以及下游區域內湍流多尺度渦結構的生成.然而對于SM 模型,圖5(a)在切片位置下游仍然能夠觀察到流向的大尺度渦結構(如圖5 中箭頭所示),結合圖4(a)可知,湍流中小尺度渦的發展可能受到近噴口過度的亞網格尺度耗散的抑制,這與文獻[7-9]對SM 模型過度耗散特點的報道一致.圖5(b)中,CKM 模型的預測結果與SM 類似,同樣可能由于該模型中局部亞網格尺度耗散較強.圖6 為展向渦 ωz的等值面,同樣采用紅色切片標記DNS 結果[22]的勢流核末端位置.可以看到,ωz的等值面在勢流段內呈平行的渦片,隨著流動向下游的發展,ωz等值面在勢流核末端附近快速卷起,隨后渦結構被拉伸并在下游發展出復雜的三維結構.圖6(a)中,SM 模型預測的渦卷起位置位于切片之后,渦片偏長,對應圖4(a)中近噴口區的不合理耗散,SM 模型抑制了流場中展向渦的發展.圖6(b),圖6(d),圖6(e) 中,CKM,LDSM 以及CSM 模型在勢流核末端位置附近都能夠觀察到展向渦的卷起(如圖6 中箭頭標記)以及下游的渦破碎,但LDSM 與CSM 模型在切片后預測的小尺度渦結構比CKM 模型豐富,再次表明CKM 模型的局部強耗散可對小尺度渦結構產生抑制.圖6(c) 中,SMSM 預測的渦卷起不顯著,但切片下游的預測結果與LDSM,CSM類似.相比SM 與CKM 模型,SMSM,LDSM 以及CSM 模型更能反映勢流核下游展向渦的多尺度特性.

圖5 對比 Q=0.1 的等值面,采用流向速度著色,紅色切片位置標記了文獻[22]中DNS 結果的勢流核末端:x=6.43Fig.5 Comparisons of the iso-surface of Q=0.1,colored by streamwize velocity.The red slice marks the end of the potential core x=6.43 of the DNS results in Ref.[22]

圖6 對比 ωz 的等值面,紅色切片位置標記了文獻[22]中DNS 結果的勢流核末端:x=6.43Fig.6 Comparisons of the iso-surface of ωz,the red slice marks the end of the potential core x=6.43 of the DNS results in Ref.[22]

2.3 速度場的POD 模態分析

對5 種SGS 模型的預測結果進行POD 分解.平行噴流主要表現為x-y平面內的二維流動特征,因此在湍流統計的流場數據中,截取三維流場的中心平面(z=1.5)為樣本,以10Δt(=0.030 75) 為采樣間隔,得到4685 個樣本.針對5 種SGS 模型,提取脈動速度分量u′(=u?〈u〉),v′(=v?〈v〉) 和w′(=w?〈w〉)的POD 模態.脈動速度u′,v′和w′的第一階模態分別對原始脈動速度場的脈動強度:〈u′u′〉,〈v′v′〉 和〈w′w′〉的貢獻最大,并在這一意義下脈動速度的第一階模態與原始脈動速度場最為接近.圖7 為流向脈動速度u′模態的累計能量貢獻率 η,可以看到,累計能量貢獻率曲線在前200 階內就達到了90%以上,這200 階模態僅為總模態數的4%,因此曲線 η 反映了POD 的快速收斂特性.但5 種SGS 模型的 η 具有不同的增長速度,SM 模型最快,而SMSM 模型最慢.脈動速度v′與w′的 η 曲線具有與圖7 相似的快速收斂特征,因此不再繪出.圖8 為脈動速度分量在前20 階模態的能量貢獻率衰減曲線 ξ.可以看到,5 種SGS 模型得出的 ξ 曲線都在前5 階內迅速衰減,說明低階模態占主導且模態表征的流場具有強的相關性.對比圖8 中第1 階模態的能量貢獻率,圖8(a)中u′模態的衰減曲線顯示SM 模型的第1 階模態貢獻率最高,SMSM 模型最低;圖8(b)中橫向速度脈動v′模態的衰減曲線顯示CKM 模型的貢獻率最高,而SMSM 模型依然排最低;圖8(c)展向速度脈動w′模態的衰減曲線中SMSM 模型的能量貢獻率排序有所提高.圖8 表明,第1 階模態的能量貢獻率已經體現出了不同模型之間的差異,接下來針對第1 階模態展開分析.

圖7 u′ 模態的累計能量貢獻率的收斂曲線Fig.7 Convergency curves of the cumulative energy contribution for the modes of u′

噴流主要表現為流向的動量傳遞.對比5 種SGS 模型預測的流向速度脈動u′的第一階POD 模態,圖9 顯示了模態幅值在中心平面(z=1.5)上的分布,采用接近零的模態幅值等值線u′=±0.003 來顯示模態的輪廓.圖9(a)中,SM 得到兩條沿流向的帶狀模態,且在噴流下游x=12～24 的范圍內一直保持高的模態幅值,對應圖5(a)中大尺度流向渦結構對流場的攪動.圖9(b)中,CKM 模型預測的模態也在下游保持沿流向的高幅值帶狀模態,在轉捩區(x=6～9)內的模態強度呈現一對低谷區,位置恰好匹配圖4(b)中的峰值耗散區位置.這表明流向速度脈動的POD 模態對亞網格尺度的黏性耗散率敏感,局部的強耗散導致脈動強度〈u′u′〉 的降低,從而改變了POD 模態的形狀.圖9(d)、圖9(e)分別為LDSM與CSM 模型預測的模態,兩者都沿流向兩側發散開來,并且在噴流下游呈現分岔或破碎的形態,模態幅值也逐漸減小,反映下游的多尺度結構以及脈動能量的逐漸消散.圖9(c)為SMSM 模型預測的模態,模態輪廓在流向上呈斷續的條狀,局部的模態幅值低,這可能與圖8(a)中SMSM 模型第一階POD 模態的能量貢獻率低有關.而在下游,SMSM 的模態呈現與CSM,LDSM 模型預測結果相似的分岔形態.因此,與SM 和CKM 模型相比,CSM,LDSM 以及SMSM 模型都合理地反映了u′模態的形態.

圖8 脈動速度分量的POD 模態的能量貢獻率Fig.8 Energy contribution rate of POD modes for the fluctuate velocity components

圖9 流向脈動速度 u′ 的第一階POD 模態云圖,實、虛等值線分別表示 u′=+0.003和u′=?0.003Fig.9 Contours of the first-order POD modes of u′,with solid contour lines for u′=+0.003 and dashed contour lines foru′=?0.003

圖10 為橫向脈動速度v′的第一階POD 模態,等值線v′=±0.003 顯示出模態的輪廓.同時,圖10中采用面內矢量 (u′,v′) 的POD 模態來顯示噴流中的流動模式.觀察圖10,v′的POD 模態呈正負相間的肋狀,并沿流向排列.在肋狀模態上下端部的附近,矢量箭頭顯示出環形的流動模式,在噴流上部的紅圓標出了逆時針環流,下部的藍圓標出了順時針環流.環流模式穿過了主流與環境流體,在圓環的下游側,高速的主流流體被送入兩側的低速環境流體中,而在環流的上游側,低速的環境流體被帶入主流流體中.同時,v′模態輪廓的橫向尺寸也沿著流向逐漸增大.可見,(u′,v′) 的POD 模態反映了噴流的流動卷吸現象[16].對比圖10(a)～ 圖10(e),SM,SMSM,LDSM 以及CSM模型預測的模態都反映了噴流的卷吸現象,其中SMSM 模型預測的旋渦尺度小,這與圖8(a) 和圖8(b) 中的低模態貢獻率對應.而圖10(b)中,CKM模型預測的模態基本沒有反映該區域內的流動卷吸,結合圖4(b),脈動速度分量u′和v′可能受到CKM模型在x=6 附近局部的強耗散抑制,導致未能預測出明顯的流動卷吸.

圖10 橫向脈動速度v'的第一階POD 模態云圖與面內矢量 (u′,v′) 的第一階POD 模態矢量圖(實、虛等值線分別表示v′=+0.003和 v′=?0.003,紅圈標示逆時針環流、藍圈標示順時針環流)Fig.10 Contours of the first-order POD modes of v' and the in plane vector (u′,v′)(Solid contour lines for v′=+0.003 and dashed contour lines for v′=?0.003.Red circle marks counterclockwise circular flow and blue circle indicates clockwise circular flow)

圖11 為展向脈動速度w′的POD 模態云圖,等值線w′=±0.003 顯示出模態的輪廓.圖11(b)～圖11(e)分別為CKM,SMSM,LDSM 以及CSM 模型預測的結果,w′的POD 模態在轉捩區(x=3～9)內呈現規則的脊狀排列,紅色表示脈動速度w′指向平面外,藍色表示脈動速度w′指向平面內,強弱相間的分布反映了w′模態對渦結構的展向拉伸,如圖11(f)所示.圖11 中的脊狀模態顯示渦拉伸主要發生在轉捩區,而在下游的模態輪廓呈碎塊狀,強度也逐漸減弱,對應渦結構逐漸發展為三維結構,展向的拉伸現象不再顯著.圖11(a)中,SM 模型預測的模態沒有反映轉捩區內的脊狀排列,在下游呈大團塊,位置與圖4(a) 中SM 模型噴口附近的過度耗散區域相對應.圖11(b)中,CKM 模型預測的模態雖然呈現了脊狀排列,但是模態輪廓在下游x=21 處仍然呈大尺度塊狀,與SM 模型的預測結果相似.相比而言,SMSM,LDSM 以及CSM 模型的預測結果則能夠合理地反映展向拉伸的流動模式.

圖11 展向脈動速度 w′ 的第一階POD 模態云圖Fig.11 Contours of the first-order POD modes of w′,with solid contour lines for w′=+0.003and dashed contour lines forw′=?0.003

圖11 展向脈動速度 w′ 的第一階POD 模態云圖 (續)Fig.11 Contours of the first-order POD modes of w′,with solid contour lines for w′=+0.003and dashed contour lines for w′=?0.003 (continued)

在計算效率方面,表1 對比了5 種SGS 模型的CPU 時間,以SM 模型為基準,LDSM,SMSM,CKM,CSM 的計算時間分別為SM 模型的:1.64,1.28,1.22 和1.00 倍,其中LDSM,SMSM 以及CKM 模型需要進行濾波運算,因此需要更多的CPU 時間,CSM 模型的計算時間最接近SM 模型.

表1 每時間步內每網格點所占用的CPU 時間Table 1 CPU time per grid point per time step

3 結論

本文開展了亞聲速可壓縮平行噴流的大渦模擬.在分析平均流動與耗散的基礎上,針對瞬時渦流以及脈動速度POD 模態表征的擬序結構,對SM,CKM,SMSM,LDSM 和CSM,5 種SGS 模型的性能進行研究,得出以下主要結論:

(1) 初步揭示了u′模態從流向的帶狀到下游分岔或破碎形態反映的多尺度特征、v′模態沿流向的肋狀排列與尺寸增長、矢量 (u′,v′) 模態的環流模式表征的流動卷吸、以及w′模態表征的流場結構在展向上受拉伸的模式.

(2) SMSM,LDSM 以及CSM 均較好地反映出湍流的多尺度特性,清晰地分辨了小尺度結構的發展過程.而SM 與CKM 模型則未能有效地預測湍流中的小尺度渦結構及部分流場特征模態.

(3) 脈動速度場的POD 模態對SGS 模型的亞網格尺度耗散敏感,表現為CKM 預測的峰值耗散區對應u′模態的低谷區因而環流模式不顯著,以及SM 未能反映轉捩區內w′模態的脊狀拉伸模式.SMSM,LDSM 和CSM 克服了以上模型的不足之處,其中CSM 模型同時兼有計算效率較高的優勢.