劉琪麟 賴煥新

(華東理工大學(xué)機(jī)械與動(dòng)力工程學(xué)院,上海 200237)

引言

擬序結(jié)構(gòu)普遍存在于湍流中[1-3],并影響湍流的發(fā)生、輸運(yùn)與耗散等過程,對(duì)湍流擬序結(jié)構(gòu)的研究、對(duì)于理解湍流機(jī)理、分析流動(dòng)相互作用有重要的意義.從大渦模擬來看,正確預(yù)測(cè)湍流擬序結(jié)構(gòu)的關(guān)鍵就在于算法和亞網(wǎng)格尺度(sub-grid-scale,SGS)模型的準(zhǔn)確性.

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

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

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

本文對(duì)亞聲速可壓縮平行噴流進(jìn)行大渦模擬.在驗(yàn)證計(jì)算結(jié)果的基礎(chǔ)上,研究SM,CKM,SMSM,LDSM 和CSM 模型的特性,分析各個(gè)模型預(yù)測(cè)的瞬時(shí)渦結(jié)構(gòu).同時(shí)采用基于SVD 的POD 分解方法提取噴流擬序結(jié)構(gòu)并分析各模型預(yù)測(cè)結(jié)果的差異,為可壓縮噴流的大渦模擬提供參考.

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

1.1 亞網(wǎng)格尺度模型

控制方程組為Favre 質(zhì)量加權(quán)濾波的三維非穩(wěn)態(tài)可壓縮N-S 方程組[21].方程中的亞網(wǎng)格尺度未知量由亞網(wǎng)格模型給出,包括亞網(wǎng)格尺度應(yīng)力τij

以及亞網(wǎng)格尺度熱通量

其中ρ,T,νsgs以及ui(i=1,2,3) 分別為密度、溫度、?；膩喚W(wǎng)格黏性系數(shù)以及速度矢量的3 個(gè)笛卡爾分量.γ 是量熱完全氣體的比熱比,空氣取γ=1.4.M為馬赫數(shù),湍流普朗特?cái)?shù)常數(shù)Prt=1[22],而,頂標(biāo)“-”表示網(wǎng)格濾波,“～”表示Favre 質(zhì)量加權(quán)的網(wǎng)格濾波,所有計(jì)算變量均為網(wǎng)格濾波后的變量,因此除非特殊說明以下省去頂部符號(hào).

Smagorinsky 模型(SM)中的亞網(wǎng)格尺度應(yīng)力為

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

選擇多尺度模型(slective mixed-scale model,SMSM)是根據(jù)局部渦矢量方向的脈動(dòng)來調(diào)整亞網(wǎng)格黏性系數(shù)

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

擬序結(jié)構(gòu)模型是利用擬序結(jié)構(gòu)與湍流中耗散的相關(guān)性來調(diào)整亞網(wǎng)格黏性系數(shù)vsgs=CΔ2|S(u)|,模型系數(shù)C=C4|FCS|3/2中常數(shù)C4=0.05[13],FCS=Q/E為擬序結(jié)構(gòu)函數(shù)，其中

為速度導(dǎo)數(shù)Dij=?ui/?xj的第二不變量,它主要由Di j中反映當(dāng)?shù)匦D(zhuǎn)的反對(duì)稱部分Aij和反映當(dāng)?shù)刈冃蔚膶?duì)稱部分Bi j組成,而為Di j的范數(shù).其中

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

擬序結(jié)構(gòu)動(dòng)能模型是在擬序結(jié)構(gòu)模型的基礎(chǔ)上,在亞網(wǎng)格黏性系數(shù)中顯式地包含亞網(wǎng)格尺度動(dòng)能ksgs,并做量綱調(diào)整后得到

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

1.2 計(jì)算對(duì)象與數(shù)值離散

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

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

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

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

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

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

1.3 POD 方法

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

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

在計(jì)算中,可對(duì)流場(chǎng)的速度分量進(jìn)行時(shí)間采樣,得到關(guān)于標(biāo)量 φ 的N個(gè)流場(chǎng)樣本,對(duì)于網(wǎng)格點(diǎn)數(shù)為M=Mx×My的二維樣本,可將每個(gè)樣本都寫成M×1的列向量,由此得到M×N的樣本矩陣

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

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

以及累計(jì)能量貢獻(xiàn)率

另外,由于各階POD 模態(tài) ψi是相互正交的,所以λ衰減越快,當(dāng)前模態(tài)所含有的能量相比下一階模態(tài)的差就越大,當(dāng)前模態(tài)所表征的物理場(chǎng)的相關(guān)性就越強(qiáng).

2 SGS 模型的特性

2.1 平均流場(chǎng)

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

圖2 LDSM 模型預(yù)測(cè)的脈動(dòng)速度無量綱功率譜Fig.2 Dimensionless power spectra of velocity signals predicted by the LDSM

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

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

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

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

2.2 瞬時(shí)流場(chǎng)的渦結(jié)構(gòu)

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

圖5 對(duì)比 Q=0.1 的等值面,采用流向速度著色,紅色切片位置標(biāo)記了文獻(xiàn)[22]中DNS 結(jié)果的勢(shì)流核末端: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 對(duì)比 ωz 的等值面,紅色切片位置標(biāo)記了文獻(xiàn)[22]中DNS 結(jié)果的勢(shì)流核末端: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 速度場(chǎng)的POD 模態(tài)分析

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

圖7 u′ 模態(tài)的累計(jì)能量貢獻(xiàn)率的收斂曲線Fig.7 Convergency curves of the cumulative energy contribution for the modes of u′

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

圖8 脈動(dòng)速度分量的POD 模態(tài)的能量貢獻(xiàn)率Fig.8 Energy contribution rate of POD modes for the fluctuate velocity components

圖9 流向脈動(dòng)速度 u′ 的第一階POD 模態(tài)云圖,實(shí)、虛等值線分別表示 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 為橫向脈動(dòng)速度v′的第一階POD 模態(tài),等值線v′=±0.003 顯示出模態(tài)的輪廓.同時(shí),圖10中采用面內(nèi)矢量 (u′,v′) 的POD 模態(tài)來顯示噴流中的流動(dòng)模式.觀察圖10,v′的POD 模態(tài)呈正負(fù)相間的肋狀,并沿流向排列.在肋狀模態(tài)上下端部的附近,矢量箭頭顯示出環(huán)形的流動(dòng)模式,在噴流上部的紅圓標(biāo)出了逆時(shí)針環(huán)流,下部的藍(lán)圓標(biāo)出了順時(shí)針環(huán)流.環(huán)流模式穿過了主流與環(huán)境流體,在圓環(huán)的下游側(cè),高速的主流流體被送入兩側(cè)的低速環(huán)境流體中,而在環(huán)流的上游側(cè),低速的環(huán)境流體被帶入主流流體中.同時(shí),v′模態(tài)輪廓的橫向尺寸也沿著流向逐漸增大.可見,(u′,v′) 的POD 模態(tài)反映了噴流的流動(dòng)卷吸現(xiàn)象[16].對(duì)比圖10(a)～ 圖10(e),SM,SMSM,LDSM 以及CSM模型預(yù)測(cè)的模態(tài)都反映了噴流的卷吸現(xiàn)象,其中SMSM 模型預(yù)測(cè)的旋渦尺度小,這與圖8(a) 和圖8(b) 中的低模態(tài)貢獻(xiàn)率對(duì)應(yīng).而圖10(b)中,CKM模型預(yù)測(cè)的模態(tài)基本沒有反映該區(qū)域內(nèi)的流動(dòng)卷吸,結(jié)合圖4(b),脈動(dòng)速度分量u′和v′可能受到CKM模型在x=6 附近局部的強(qiáng)耗散抑制,導(dǎo)致未能預(yù)測(cè)出明顯的流動(dòng)卷吸.

圖10 橫向脈動(dòng)速度v'的第一階POD 模態(tài)云圖與面內(nèi)矢量 (u′,v′) 的第一階POD 模態(tài)矢量圖(實(shí)、虛等值線分別表示v′=+0.003和 v′=?0.003,紅圈標(biāo)示逆時(shí)針環(huán)流、藍(lán)圈標(biāo)示順時(shí)針環(huán)流)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 為展向脈動(dòng)速度w′的POD 模態(tài)云圖,等值線w′=±0.003 顯示出模態(tài)的輪廓.圖11(b)～圖11(e)分別為CKM,SMSM,LDSM 以及CSM 模型預(yù)測(cè)的結(jié)果,w′的POD 模態(tài)在轉(zhuǎn)捩區(qū)(x=3～9)內(nèi)呈現(xiàn)規(guī)則的脊?fàn)钆帕?紅色表示脈動(dòng)速度w′指向平面外,藍(lán)色表示脈動(dòng)速度w′指向平面內(nèi),強(qiáng)弱相間的分布反映了w′模態(tài)對(duì)渦結(jié)構(gòu)的展向拉伸,如圖11(f)所示.圖11 中的脊?fàn)钅B(tài)顯示渦拉伸主要發(fā)生在轉(zhuǎn)捩區(qū),而在下游的模態(tài)輪廓呈碎塊狀,強(qiáng)度也逐漸減弱,對(duì)應(yīng)渦結(jié)構(gòu)逐漸發(fā)展為三維結(jié)構(gòu),展向的拉伸現(xiàn)象不再顯著.圖11(a)中,SM 模型預(yù)測(cè)的模態(tài)沒有反映轉(zhuǎn)捩區(qū)內(nèi)的脊?fàn)钆帕?在下游呈大團(tuán)塊,位置與圖4(a) 中SM 模型噴口附近的過度耗散區(qū)域相對(duì)應(yīng).圖11(b)中,CKM 模型預(yù)測(cè)的模態(tài)雖然呈現(xiàn)了脊?fàn)钆帕?但是模態(tài)輪廓在下游x=21 處仍然呈大尺度塊狀,與SM 模型的預(yù)測(cè)結(jié)果相似.相比而言,SMSM,LDSM 以及CSM 模型的預(yù)測(cè)結(jié)果則能夠合理地反映展向拉伸的流動(dòng)模式.

圖11 展向脈動(dòng)速度 w′ 的第一階POD 模態(tài)云圖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 展向脈動(dòng)速度 w′ 的第一階POD 模態(tài)云圖 (續(xù))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)

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

表1 每時(shí)間步內(nèi)每網(wǎng)格點(diǎn)所占用的CPU 時(shí)間Table 1 CPU time per grid point per time step

3 結(jié)論

本文開展了亞聲速可壓縮平行噴流的大渦模擬.在分析平均流動(dòng)與耗散的基礎(chǔ)上,針對(duì)瞬時(shí)渦流以及脈動(dòng)速度POD 模態(tài)表征的擬序結(jié)構(gòu),對(duì)SM,CKM,SMSM,LDSM 和CSM,5 種SGS 模型的性能進(jìn)行研究,得出以下主要結(jié)論:

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

(2) SMSM,LDSM 以及CSM 均較好地反映出湍流的多尺度特性,清晰地分辨了小尺度結(jié)構(gòu)的發(fā)展過程.而SM 與CKM 模型則未能有效地預(yù)測(cè)湍流中的小尺度渦結(jié)構(gòu)及部分流場(chǎng)特征模態(tài).

(3) 脈動(dòng)速度場(chǎng)的POD 模態(tài)對(duì)SGS 模型的亞網(wǎng)格尺度耗散敏感,表現(xiàn)為CKM 預(yù)測(cè)的峰值耗散區(qū)對(duì)應(yīng)u′模態(tài)的低谷區(qū)因而環(huán)流模式不顯著,以及SM 未能反映轉(zhuǎn)捩區(qū)內(nèi)w′模態(tài)的脊?fàn)罾炷Ｊ?SMSM,LDSM 和CSM 克服了以上模型的不足之處,其中CSM 模型同時(shí)兼有計(jì)算效率較高的優(yōu)勢(shì).