三維效應對超聲速湍流燃燒流動大渦模擬的影響研究

劉宏鵬

國防科技大學 空天科學學院，湖南 長沙410073

近幾十年來，隨著電子計算機在運算速度與數據存儲技術上的迅速發展，以及數值計算方法的進步，計算流體力學（CFD）已被廣泛應用到超聲速湍流燃燒流動的基礎研究和工程設計中[1-3]。在數值模擬研究中，當前普遍采用的是雷諾平均（RANS）方法，該方法早已成功用于工程實際的復雜流動模擬，但只能給出流場的統計平均量，無法給出流動細節。直接數值模擬（DNS）雖然可給出全時、空域的詳細瞬時流動過程，但由于受到計算機條件與計算方法的制約，在現階段只適用于具有簡單幾何邊界的較低或中等雷諾數下的湍流流動模擬。近些年來，介于RANS和DNS之間的大渦模擬（LES）方法廣被應用于超聲速湍流燃燒流動的模擬研究中[4-14]。

為節約計算量，國內外均出現應用二維LES研究超聲速湍流燃燒的工作。如D.Chakraborty等[15]采用二維隱式LES 研究了不同化學反應機理對超聲速湍流燃燒混合層數值模擬結果的影響。Y.L.Zhang等[16]、C.Qian等[17]分別對超聲速湍流燃燒混合層流動開展了二維LES研究，重點觀察了火焰的點火、傳播、熄火及火焰不穩定等現象。Z. X.Ren等[18-19]采用二維LES研究了超聲速湍流燃燒多項流中漩渦、激波和燃燒放熱等物理化學現象的相互干擾。這些研究可定性地揭示部分湍流的物理特性，但在實際應用中難以提出可靠的數據。三維屬性可對湍流場的數值模擬產生重要影響：P.A.Davidson 在其論著[20]中指出，三維湍流中至關重要的漩渦拉伸作用在二維流動中將失效，由此導致二維湍流呈現出獨特特性；此外，M.Breuer[21]具體對比了繞圓柱不可壓縮流動的二維/三維LES結果，發現忽略三維效應后，圓柱尾流區的速度分布、分離區大小以及圓柱物面的壓強（壓力）、摩阻系數分布以及背壓系數均產生較大誤差。盡管三維效應對于湍流非定常數值模擬極其重要已成為共識，但以往研究大多集中于分析三維效應對于無反應湍流流動數值模擬準確性的影響及其機理。在超聲速湍流燃燒領域，由于高馬赫數湍流和燃燒間的強烈耦合相互作用，二維LES模擬導致的湍流場結果的誤差必然影響燃燒場的數值模擬結果的準確性，使得二維湍流燃燒場呈現獨特特性，而關于二維/三維超聲速湍流燃燒LES 模擬結果的詳細對比研究并不多見。

綜上所述，本文旨在探究三維效應對于超聲速湍流燃燒流動LES 數值模擬準確性的影響。本文介紹了控制方程，對一個對流馬赫數Mac=0.61 的超聲速湍流燃燒混合層流動分別開展了二維/三維LES研究，通過對比二維/三維結果分析了三維效應對LES準確模擬超聲速湍流燃燒流動的重要性。

1 計算方法

笛卡兒直角坐標系下的考慮多組分燃燒化學反應的質量加權濾波后的Navier-Stokes（N-S）控制方程組包含組分輸運、質量守恒、動量守恒和能量守恒方程，分別為

式中：ρ，ui，p與T分別為氣體的密度、速度、壓力與靜溫；E為單位質量總能量

式中：ns為組分總數，Ms、Ys、Ds、hs與ω?s分別為組分s的摩爾質量、質量分數、質量擴散系數、單位質量的絕對焓值以及單位時間、單位體積的化學反應質量源項，式(1)中的Sct,s為組分s的湍流施密特數，在本研究中取為0.9。需要進一步說明的是，為確保組分質量分數之和為1，通過1 減去其余組分質量分數之和而求得最后一組分的質量分數，即

在本文中，采用氫氣-空氣的7 組分(H2, O2, H2O, H,O, OH, N2) 8 反應模型[22]化學反應機理模擬燃料為氫氣的超聲速湍流燃燒混合層流動，最后一組分選定為氮氣N2。?分別為分子黏性應力項和熱量擴散項

式中：μ和μt分別為為混合氣體分子[動力]黏度和亞格子渦黏度，λ為混合氣體的熱傳導系數，Prt為湍流普朗特數，取為0.5。此外，通過分子動理論[23]獲得各組分的分子黏度、熱傳導系數、質量擴散系數，進而通過Wilke的混合律[24]得到混合氣體的分子黏度、熱傳導系數、質量擴散系數。采用關于溫度的多項式擬合獲得各組分的定壓比熱容和靜焓，擬合系數見參考文獻[25]。為封閉控制方程組，需要加入氣體狀態方程

式中：M為混合氣體的相對分子質量，Ru=8314J/(kmol·K)為普適氣體常數。為亞格子雷諾應力

通過Smagorinsky模型模化亞格子雷諾應力，公式為

2 超聲速湍流燃燒混合層流動二維/三維LES對比研究

2.1 算例描述及其驗證

超聲速湍流燃燒混合層流動的來流條件[16]選取典型的超燃沖壓發動機燃燒室內來流條件，見表1，其中空氣來流馬赫數為3.0，溫度為1100K，燃料選為氫氣和氮氣的混合物，以接近聲速的速度（Ma=1.13）與空氣摻混。混合層流動示意圖如圖1 所示，其中x，y，z分別表示流向、法向和展向。對流馬赫數定義為其中下標1和2分別表示空氣和燃料來流，c表示聲速。通過表1計算可知本文研究的混合層Mac=0.61。此外，基于入口動量厚度的雷諾數Reθ=2226。

表1 超聲速混合層燃燒數值試驗自由來流條件Table1 Freestream conditions for computational experiment of supersonic turbulent combustion for a mixing layer flow

圖1 超聲速混合層燃燒數值試驗流動示意圖Fig.1 Schematic of supersonic turbulent combustionflow for a mixing layer flow

入口邊界配置雙曲正切函數型[26]的流向速度分布：

式中：b為控制參數，取為4×10-4m，法向v和展向速度w均為0。為刺激失穩，在入口邊界處除了配置時均速度外，額外增加白噪聲擾動[26]的脈動速度分量，擾動幅值為0.01uair。此外，如圖1 所示，出口為超聲速出口邊界條件，展向采用周期邊界條件。

二維LES 所采用的網格計算域尺寸為1.2m（流向）×0.3m（法向），網格分布為1024（流向）×256（法向），流向網格均勻分布，法向網格在空氣流和氫氣流的交接處加密。將二維網格向展向拉伸0.15m，并均布128個網格點獲得三維網格。網格示意圖如圖2所示。

圖2 x-y平面網格示意圖(為使圖清晰，x、y方向網格點數均縮小一倍)Fig.2 Schematic of computational mesh at x-y center plane(Every 2th point is displayed at both streamwise and transverse directions for clarity)

對流無黏格式采用5 階PPWEON-Z 格式[27]，黏性項采用2階中心差分格式進行離散。時間推進采用雙時間隱式迭代方法[28]，時間步長設定為1×10-7s，內迭代步數為5 步。數值模擬采用氫氣-空氣的7 組分(H2, O2, H2O, H, O, OH,N2)8反應模型[22]。首先計算推進兩倍最大通流時間T=2Lx/(U1+U2)得到充分發展湍流流場，然后統計6 倍最大通流時間得到湍流統計相關量[29]，其中Lx=1.2m 為流向計算域尺度。作者在參考文獻[27]將三維LES結果與DNS結果進行了對比，發現LES計算獲得的雷諾切應力及流向、法向脈動速度均方根分布等與DNS結果相當，驗證了本文所采用的數值方法的可靠性和準確性。

2.2 流動和火焰結構

基于氫原子H的混合分數定義為

化學當量值Zst=1/(1+S)=0.088,其中S=rYH2,F/YO2,A=10.345。在該算例中，YH2,F=0.3 為燃料一側來流中氫氣的質量分數，YO2,A=0.232為空氣一側來流中氧氣的質量分數，r=8為H2/O2充分燃燒時的質量比。圖3對比了二維/三維LES的瞬時溫度分布，其中溫度云圖中的黑色曲線為ZH=Zst等值線。由圖3可見，無論二維還是三維LES，燃燒均主要發生在ZH=Zst周圍。然而，二維LES 下的流動結構以及火焰形態卻明顯與三維LES結果不同：在流動結構上，二維LES計算獲得的流動圖像呈現明顯的旋渦卷起、分裂和破碎等變化，與Zhang 等[16]、Wang等[30]、Lele等[31]研究采用二維LES/DNS獲得的流動圖像一致，而三維LES計算獲得的湍流流動結構呈不規則狀，與Edwards等[32]的三維RANS/LES結果相似；在火焰形態上，二維LES的火焰偏薄，呈條帶狀，而三維LES 的火焰呈厚的擴散狀。圖4對比了二維/三維LES 瞬時H2O 質量分數分布。由圖4 可見，二維/三維LES 的自點火位置大致相同，約位于x=0.5m 處，且瞬時H2O的質量分數峰值均在0.2以上。

圖3 二維/三維LES的瞬時溫度分布Fig.3 Instantaneous temperature contours for 2D/3D LES

圖4 二維/三維LES的瞬時H2O質量分數分布Fig.4 Instantaneous mass fraction of H2O contours for 2D/3D LES

2.3 湍流旋渦與火焰相互干擾

圖5 對比了二維/三維LES 下的渦量分布。由圖5 可見，二維/三維LES 的渦量形態顯著不同：二維LES 的渦量呈明顯的螺旋狀表征渦的生成、卷起、并對和破碎；而三維LES渦量呈不規則狀。可壓縮情況下（忽略外力）的渦量方程的張量形式為

圖5 二維/三維LES的渦量Fig.5 Vortex contours for 2D/3D LES

式中：DΩ/Dt表示隨體導數，Ω=▽×u為渦量，μv為第二黏度。為簡化討論，假設流體不可壓且在二維情況下，式(16)退化為

式(17)表明此時渦量的演化僅受擴散作用的影響。渦量演變的重要支配項，即渦管的拉伸和壓縮機制(Ω??)u由于空間維度自三維退化為二維而失效，因而造成了圖5所示的二維/三維湍流渦量場的顯著區別。

圖6 顯示了二維/三維LES 計算獲得的渦量、質量分數及OH 化學反應速率的等值線。其中，質量分數等值線范圍為0.008 ≤ZH≤0.4，包含化學當量值Zst=0.088，相應也包含高溫火焰區域。OH 生成速率等值線范圍取為ω?OH=10～100(kg/m3)/s，為高OH 生成速率區域。由圖6 可見，二維/三維情況下，高溫火焰區域與高OH生成速率區域一致。火焰形態受渦量場影響而呈現不同形態：二維情況下，受旋渦運動影響，燃燒場也呈現出渦的卷起、并對和破碎等現象；而三維情況下，燃燒場呈不規則狀。

圖6 二維/三維LES的湍流旋渦和火焰相互干擾Fig.6 Flow pattern for interactions between turbulent vortex and combustion flames

2.4 混合層厚度增長率

定義δ0、δω和δθ分別為邊界層名義厚度、渦量厚度和動量厚度，表達式如下

式中：y_up 和y_down分別表示速度U1-0.1ΔU(無量綱值為0.9)和U2+0.1ΔU(無量綱值為0.1)所對應的y向坐標。在動量厚度的定義中，積分上下限通常分別為Ly/2 和-Ly/2。而在當前研究中，其分別以y_up 和y_down替代，原因是空間發展的超聲速湍流燃燒混合層流動的點火將誘導激波，激波使得混合層外流動速度降低而與來流條件不一致，混合層分布[y_up,y_down]以外的動量厚度積分將“污染”整體積分值。

圖7 給出了二維/三維LES 計算獲得的混合層厚度沿流向發展結果，各厚度均以入口處的值作無量綱化。由圖7可見，二維/三維LES計算得到的三種混合層厚度具有相同的發展趨勢，在點火位置（約位于x/δω0=375）前有一段近似線性的發展區域，點火位置附近區域混合層厚度呈非線性發展，在點火位置后同樣有一段近似線性的發展區。點火位置后的混合層厚度近似線性發展區域中，無量綱湍動能分布在不同流向站位處的趨于一致，如圖8所示。因此，可判定該區域為流動自相似區。采用最小二乘法對流動自相似區的混合層厚度進行線性擬合，兩種模式給出的三種混合層厚度增長率見表2。由表2 可知，二維LES 給出的自相似區三種混合層厚度增長率均明顯低于三維LES結果。

圖7 二維/三維LES的混合層厚度沿流向發展Fig.7 Thickness development of mixing layer flows for 2D/3D LES

圖8 超聲速湍流燃燒混合層流動自相似區判定Fig.8 Determination of self-similar regions for the mixing layer flows

表2 二維/三維大渦模擬下混合層增長率Table 2 Mixing layer thickness growth rate for 2D/3D LES

2.5 速度型、雷諾應力和湍動能

對自相似段幾個不同流向站位分布做算術平均后，圖9給出了二維/三維LES和DNS[33-34]的自相似段湍流統計相關量的分布特征。由圖9可見，三維LES的速度型、湍動能及雷諾剪切應力均與DNS 結果相當。二維/三維LES 的無量綱流向速度分布基本相符。然而，相較于三維LES結果，二維LES 的湍動能峰值增大了約59%，而雷諾剪切應力峰值減小了約62%。此外，二維LES 的湍動能和雷諾剪切應力在法向的擴散更加強烈，導致其分布更加寬闊。Liou等[35]的無反應可壓縮二維LES下的混合層流動的湍動能和雷諾剪切應力同樣呈現當前結果的分布特征。

圖9 二維/三維LES的統計平均速度型、湍動能和雷諾剪切應力對比Fig.9 Mean profiles of streamwise velocity,turbulent kinetic energy and Reynolds shear stress for 2D/3D LES compared to DNS results

2.6 湍流能譜

圖10對比了二維/三維LES的能譜。其中，斯特勞哈爾（Strouhal）數定義為：Sr=f·δθ0/Uc，δθ0表示入口邊界層厚度，而Uc= (U1+U2)/2，E(k) = 0.5[E(u) + E(v) + E(w)]。在圖10（b）中，為清晰起見，將E(v2)除去了因子102。觀測點位于混合層中心位置x=0.83Lx，y*=(y-yc)/δω=-0.46（三維），-0.17（二維）。由圖10可見，二維和三維情況下的能譜均具備寬闊的頻域范圍，表明流動達到充分發展湍流狀態。此外，二維/三維LES 計算獲得的在含能和耗散尺度之間的區間能譜斜率均滿足-5/3律。然而，由圖10（b）的法向速度v的能譜可知，二維LES 計算獲得的法向速度在各個尺度上的脈動相較于三維LES結果都更加強烈。

圖10 二維/三維LES的能譜對比Fig.10 Energy spectrum for 2D/3D LES

2.7 統計平均溫度及其脈動

圖11 對比了二維/三維LES 所給出的統計平均溫度云圖和計算域流向出口溫度分布。由圖11（a）、圖11（b）可見，三維LES計算所得火焰溫度明顯高于二維LES結果。截取計算域出口溫度分布得到圖11（c）：三維LES計算域出口處溫度峰值達到2.74T0(2040K)，而二維LES 的峰值僅為1.99T0(1480K)，比前者減小了27%。以上結果表明，二維LES所預測的統計平均溫度被嚴重低估。

圖12給出了二維/三維LES計算得到的溫度脈動結果，其中圖12（c）給出了溫度脈動峰值及峰值所在法向位置沿流向的分布。由圖12（a）、圖12（b）可見，二維計算得到的溫度脈動均方根峰值明顯高于其對應的三維結果。二維/三維LES 均呈現點火導致最大溫度脈動，如圖12（c）所示的結果。但點火后，二維/三維LES的所取得的峰值呈現截然相反的對稱性：前者趨于燃料一側（即yTrms/δω0>0），而后者趨于氧化劑一側（即yTrms/δω0＜0）。

圖12 二/三維LES的溫度脈動結果Fig.12 Statistical temperature fluctuation for 2D/3D LES

正如2.3節所揭示的，二維LES的湍流旋渦運動中的渦拉伸、壓縮機制失效，導致二維LES的旋渦結構與三維LES顯著不同，燃料和氧化劑相遇的火焰前鋒接觸面面積從而存在差別，進而致使二維LES難以獲得準確的燃料消耗率/燃燒效率，最終導致圖11所示的二維LES的統計平均火焰溫度顯著低于三維LES結果。

圖11 二維/三維LES的統計平均溫度Fig.11 Mean temperature for 2D/3D LES

3 結論

本文通過對對流馬赫數為0.61的超聲速湍流燃燒混合層流動開展二維/三維大渦模擬，探討了三維效應對于超聲速湍流燃燒流動LES的重要性。研究表明：

（1）相較于三維LES，二維LES下的超聲速湍流燃燒混合層的增長率大幅減小，其中動量厚度減小幅度約為23%。二維LES的無量綱速度型與三維LES結果一致，雷諾剪切應力明顯偏小，峰值減小約62%；湍動能明顯偏大，峰值增大約59%。此外，雷諾剪切應力和湍動能的擴散明顯增強，導致其分布寬度增加。此外，二維LES 計算所得的湍動能能譜與三維LES 結果相似，在含能渦區域均滿足-5/3分布律。

（2）由于二維流動中漩渦拉伸、壓縮機制失效，致使二維LES計算所得的湍流流動結構與三維LES存在顯著差別：二維LES中湍流呈現清晰的旋渦生成、發展、并對和破碎演化規律，而三維LES的湍流結構呈不規則狀。由于湍流和燃燒的相互作用，二維/三維的燃燒場（包括混合分數場、溫度場等）也存在顯著差別。旋渦運動對燃料和氧化劑進行摻混，摻混界面為火焰前鋒面。由于二維/三維流場的旋渦結構不同，致使火焰前鋒面面積不同，導致二維LES 計算獲得的統計平均火焰溫度顯著低于三維LES 結果。在本算例中，計算域出口處的二維LES 的火焰溫度峰值比三維LES 結果低約27%。此外，二維LES的溫度脈動顯著高于三維LES結果。

綜上所述，三維效應對于LES 準確模擬超聲速湍流燃燒流動至關重要。