可壓縮槽道流動的大渦模擬研究

王鎖柱，楊天鵬，楊依峰，秦緒國，蘇 偉

可壓縮槽道流動的大渦模擬研究

王鎖柱，楊天鵬，楊依峰，秦緒國，蘇 偉

（北京航天長征飛行器研究所，北京，100076）

可壓縮湍流邊界層是高速飛行器研究面臨的基礎科學問題，對其流動機理的認識對于高速飛行器氣動力熱設計具有重要意義。采用大渦模擬方法對可壓縮槽道流動進行了數值模擬研究，并基于Favré過濾建立了適用于可壓縮流動的大渦模擬方法，將混合亞格子模型延拓到了動態形式，形成了適用于可壓縮流動的動態混合模型。通過分析亞格子模型、網格疏密和計算格式對數值模擬結果的影響，結果表明所建立的大渦模擬方法在較少網格數下能夠獲得可靠的計算結果。數值模擬還獲得了可壓縮槽道流動時間轉捩的過程，捕捉到了大尺度旋渦結構和近壁速度條帶結構的演化過程。

可壓縮湍流；槽道流動；大渦模擬

0 引 言

可壓縮湍流現象廣泛存在于高速飛行器的邊界層流動中，具有多尺度、高頻脈動、強非線性等特點，對其流動機理的研究對于高速飛行器設計具有重要意義。隨著計算流體力學（Computational Fluid Dynamics, CFD）和計算機技術的迅速發展，數值模擬方法已經成為可壓縮湍流問題的重要研究手段。

大渦模擬（Large Eddy Simulation，LES）是綜合考慮了工程應用的要求以及計算能力的局限而發展起來的一種介于雷諾平均方法（Reynolds Averaged Navier-Stokes，RANS）和直接數值模擬（Direct Numerical Simulation，DNS）之間的數值方法。一方面，LES比RANS易于建立更為普適通用的湍流模型，能夠獲得真實的瞬態流場信息，可用于預測分離、轉捩等復雜流動。另一方面，有效的亞格子模型可使LES所需網格總數比DNS少一到兩個數量級，從而其計算量要明顯少于DNS。亞格子模型是大渦模擬的核心問題之一，對于大渦模擬的準確性和計算效率起關鍵作用。Smagorinsky模型[1]是最早的亞格子模型，其簡單易行，計算穩定性和魯棒性也較好，但包含非普適的經驗常數，同時耗散性過大。隨后發展的尺度相似模型[2]能夠較準確地表達大尺度和小尺度間的動量輸運關系，其缺陷是耗散嚴重不足，數值計算存在穩定性問題。綜合Smagorinsky模型和尺度相似模型的優點可形成混合模型[3]，既有正確的亞格子動量輸運，也有足夠的亞格子耗散。Germano等[4]提出了著名的動態模型，該模型本身不提出新的模型，而是建立在給定的基準模型上，通過動態方法確定模型中的系數，如動態Smagorinsky模型、動態混合模型[5]等。

可壓縮槽道流動作為一種可壓縮壁湍流，不僅包含了壁面效應，還可體現可壓縮效應的影響，是一種典型的可壓縮湍流流動。Gamet等[6]對馬赫數為0.2的槽道流動進行了DNS研究，由于馬赫數較低，流場特征與不可壓湍流非常接近。Lenormand等[7]采用LES對馬赫數為0.5和1.5的槽道湍流進行了模擬，研究了不同亞格子模型對計算結果的影響。Mossi和Sagaut[8]也對馬赫數為0.5和1.5的槽道湍流進行了模擬，分析了不同激波捕獲格式對LES結果的影響。李新亮等[9]采用DNS研究了馬赫數為0.8的可壓縮槽道湍流場，分析了壓縮性效應對近壁相干結構的影響。

本文采用大渦模擬方法對馬赫數為0.5的可壓縮槽道流動進行了數值模擬研究，分析了亞格子模型、網格疏密和計算格式對大渦模擬結果的影響，并給出了可壓縮槽道流動時間轉捩過程中大尺度流向渦和近壁速度條帶結構的演化過程。

1 物理模型和數值方法

1.1 物理模型

槽道流動示意如圖1所示，兩個平行平板間充滿粘性可壓縮流體沿流向運動。設笛卡爾坐標系的坐標軸、和方向分別平行于流向、法向和展向，為槽道半高。本文取槽道的半高=1，=1對應槽道的中心線。

圖1 槽道流動示意

在不可壓縮槽道流動的數值模擬中，為了保證流向的周期性，通常假設流動是由均勻壓力梯度驅動的，即以流向壓差平衡上下壁面的摩擦。但在可壓縮槽道流動的模擬中，若仍假設由平均壓力梯度驅動，則流場中必然誘導出密度和溫度梯度，進而影響整個流場參變量的分布，導致流向的周期條件無法維持。本文采用Coleman等[10]提出的可壓縮槽道流動的驅動模型，認為流體是在均勻體積力的作用下流動，即通過在流向動量方程中添加體積力項1對壁面摩擦力進行平衡，以保證流場變量滿足流向的周期性條件，同時在能量方程也需要引入一項11。為避免計算體積力時出現數值不穩定，本文采用Lenormand等[11]發展的可壓縮流動修正算法來更新每一時間步的體積力。

1.2 大渦模擬方法

引入Favré過濾后，除密度和壓強仍采用直接空間過濾表示外，其余流場變量，如速度場、溫度等均取Favré過濾量，則可導出基于Favré過濾后的N-S方程組為

其中，模型系數、I可基于瞬態流場動態計算得到。

采用有限差分法對大渦模擬控制方程進行離散求解。對流項經通量分裂后采用五階迎風型緊致格式[12]，粘性項則采用六階中心型緊致格式[13]，時間離散采用滿足TVD特性的三階Runge-Kutta方法。

1.3 初始條件和邊界條件

對于充分發展的槽道湍流，流向和展向是均勻同性的。因此可以沿流向和展向采用周期邊界條件，法向上下壁面采用等溫無滑移條件，即=w和w=0。

槽道初始的速度場、溫度場按二維層流Poiseuille流動的解析解給出，密度場則設為均勻分布。為了誘導初始層流流動順利轉捩為充分發展的湍流流動，在流向速度中需要加入一定幅值的隨機擾動。具體初始條件如下：

2 結果分析

本文對馬赫數為0.5的亞聲速槽道流動的時間轉捩過程進行了大渦模擬研究。計算域沿流向、法向和展向的長度分別為2π、2和4π/3，基于體積平均密度、體積速度、槽道半高和壁面粘性系數的雷諾數為3000。計算網格在流向和展向均勻分布，在法向采用雙曲拉伸對壁面網格加密。槽道流動在流向隨機擾動的作用下發生了明顯的時間轉捩，從層流轉變為完全湍流。本文對處于完全湍流狀態下的流場進行統計平均，由于流向和展向的均勻性，統計除了按時間平均，也沿流向和展向進行空間平均。

2.1 亞格子模型影響

為了考察不同亞格子模型對計算結果的影響，本文采用4種不同亞格子模型分別對可壓縮槽道流動進行了數值模擬。表1給出了不同亞格子模型計算得到的典型平均量的對比。表中f為壁面摩擦系數、τ為摩擦雷諾數、u為壁面摩擦速度、c為槽道中心處速度、b為體積速度。DMM代表本文所發展的動態混合模型，SM代表Smagorinsky模型，DSM代表動態Smagorinsky模型，HSSM代表混合模型。表中還給出了幾組對比數據，NAH代表Niederschulte等[14]的實驗測量結果，KMM代表Kim等[15]的不可壓縮DNS數據，Dean為經驗公式給出的結果[16]。可以發現，SM和DSM計算得到的壁面摩擦系數和壁面摩擦速度相對于DMM誤差稍大，而HSSM的結果與KMM很接近。HSSM計算得到的槽道中心速度較為準確，其他3種模型的結果都稍大。

圖2給出了采用不同亞格子模型計算得到的流向平均速度分布。在全局坐標軸下，4種模型得到的速度型相差不大，SM和DSM的速度型在靠近槽道中心時偏高，而HSSM的結果則偏低。在壁面坐標軸+下，4種模型在粘性底層得到的結果都較為一致，而在對數區SM和DSM的速度分布出現了上抬，因為該兩種模型所計算得到的壁面摩擦速度稍小。

表1 不同亞格子模型計算的平均量對比

Tab.1 Mean Flow Variables for Different SGS Models

算例CfReτuτucuc/uτub/uτ NAH8.29×10-31790.06441.1517.8615.54 KMM8.18×10-31800.06371.1618.2015.63 Dean8.29×10-3——1.16—— DMM8.32×10-31960.06371.1818.5215.92 SM7.96×10-31910.06241.1818.9116.28 DSM7.76×10-31880.06181.1819.0916.53 HSSM8.19×10-31950.06281.1618.4715.94

圖2 不同亞格子模型計算的流向平均速度

圖3給出了采用不同亞格子模型計算得到的速度脈動均方根分布。由圖3可以看出，SM和DSM得到的流向速度脈動均方根的峰值要比DMM和KMM的結果大，峰值位置也有所偏離；而HSSM的峰值偏小。對于法向速度脈動均方根，相對于DMM的結果，SM、DSM和HSSM得到的峰值依次減小，峰值位置也向槽道中心偏離。展向速度脈動均方根的分布與法向類似。可以看出，本文所建立的動態混合模型在總體上要優于其他3種模型，得到的數值模擬結果也更為準確。

圖3 不同亞格子模型計算的速度脈動均方根

2.2 網格疏密影響

為了考察網格疏密對計算結果的影響，對網格加密后進行數值模擬，亞格子模型均采用動態混合模型。表2給出了網格特征參數對比，n、n、n分別表示流向、法向和展向的網格點數，Δx和Δz是壁面坐標系下的流向、展向網格間距，min(Δy)是壁面坐標系下法向離壁面第1層的網格間距，max(Δy)則是槽道中心的網格間距。=0.5為原始網格，=0.5為對法向網格進行加密，0.5為對展向網格進行加密。

表3給出了網格加密后計算得到的典型平均量對比。由表3可以看出，法向網格加密后得到的壁面摩擦系數和壁面摩擦速度略小于加密前的量，而展向網格加密后的結果則變大。法向或展向加密計算得到的槽道中心速度均為1.15，略小于加密前的結果。

表2 網格參數對比

Tab.2 Comparison of Grid Characteristics

算例nxnynzΔx+min(Δy+)max(Δy+)Δz+ KMM192160129120.054.47 Ma=0.541656530.751.0512.312.81 Ma=0.5-y411296530.770.516.212.82 Ma=0.5-z416512931.521.0812.66.57

表3 不同網格計算的平均量對比

Tab.3 Mean Flow Variables for Different Grids

算例CfReτuτucuc/uτub/uτ NAH8.29×10-31790.06441.1517.8615.54 KMM8.18×10-31800.06371.1618.2015.63 Dean8.29×10-3——1.16—— Ma=0.58.32×10-31960.06371.1818.5215.92 Ma=0.5-y8.18×10-31960.06271.1518.3415.91 Ma=0.5-z8.59×10-32000.06421.1517.9115.53

圖4給出了網格加密后計算得到的流向平均速度分布。在全局坐標系下，網格加密后的速度型相差不大，只是在靠近槽道中心處偏低，與KMM的誤差也很小。在壁面坐標系下，法向加密后的分布與加密前的結果較為一致，而展向加密后的速度型在對數區略有降低。

圖4 不同網格計算的流向平均速度

圖5給出了網格加密后計算得到的速度脈動均方根分布。由圖5可以看出，法向加密后對流向速度脈動均方根的分布影響不大，而展向加密后分布的峰值有所減小，但峰值位置無變化。對于法向速度脈動均方根，法向加密后得到的分布基本不變，展向加密后在近壁面處與KMM符合得更好，而在靠近槽道中心時誤差稍大。網格加密后展向速度脈動均方根的分布與法向相似。根據以上分析，綜合考慮數值模擬精度和計算量的影響，本文采用的網格參數已可滿足槽道流動大渦模擬的需求。

圖5 不同網格計算的速度脈動均方根

2.3 計算格式影響

為了考察普通迎風型格式是否可用于進行大渦模擬計算，本文嘗試采用目前較流行的幾種迎風型格式開展了數值模擬研究，包括Harten-Yee TVD格式、AUSMDV格式和Roe格式。但在同樣的網格參數下，即使將流向速度擾動幅值增大到0.5，流動也未能發生轉捩，始終保持層流狀態。顯然，以上幾種迎風型格式的數值耗散過大，不適用于湍流機理的大渦模擬研究。

2.4 流場結構分析

大渦模擬的一大優點在于能夠獲得大量瞬態流場，可以從中提取大尺度相干結構隨時間或空間的演化過程，有利于對流動現象的深化認識。速度梯度張量的第二不變量是流向渦和展向渦的綜合體現，本文采用法則來識別旋渦結構，即主要觀察速度梯度張量第二不變量的等值面。

圖6給出了計算過程中3個典型時刻=0.85的等值面分布。由圖6可以看出，在計算的初始階段（=5時）流場中幾乎不存在大尺度相干結構；=50時，流場結構非常紊亂，存在大量不規則小渦結構；隨著時間的進一步推移，開始出現大尺度旋渦結構，到=100時流場中已存在大量的準流向渦，這些準流向渦與壁面和流向均呈一定角度，并在向下游運動的過程中不斷演化發展。

圖6 不同時刻大尺度相干結構

圖7給出了對應時刻瞬態流場=0.04法向截面的流向速度云圖。由圖7可知=5時流向速度較小，分布也較均勻；當=50時，流向速度出現大幅值的擾動，但未出現明顯的展向速度條帶結構；=100時，沿展向開始出現明顯的高速區和低速區，并且呈條帶狀沿展向交錯分布。這樣的近壁速度條帶結構是壁湍流重要的擬序結構。

圖7 不同時刻法向截面流向速度云圖

3 結束語

本文采用大渦模擬方法對可壓縮槽道流動進行了數值模擬研究，分析了亞格子模型、網格疏密以及計算格式對數值模擬結果的影響。結果表明，本文采用的大渦模擬方法和發展的動態混合模型能夠適用于可壓縮湍流的數值模擬，采用較少的網格可以獲得準確的結果。數值模擬還獲得了可壓縮槽道流動時間轉捩的過程，捕捉到了大尺度旋渦結構和近壁速度條帶結構的演化過程。

[1] Smagorinsky J. General circulation experiment with the primitive equations[J]. Monthly Weather Review, 1963, 91(3): 99-165.

[2] Bardina J, Ferziger J H, Reynolds W C. Improved subgrid scale models for large-eddy simulation[R]. AIAA Paper 80-1357, 1980.

[3] Meneveau C, Katz J. Scale-invariance and turbulence models for large eddy simulation[J]. Annual Review of Fluid Mechanics, 2000(32): 1-32.

[4] Germano M, et al. A dynamic subgrid-scale eddy viscosity model[J]. Physics of Fluids A, 1991, 3(7): 1760-1765.

[5] Zang Y, Street R L, Koseff J R. A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows[J]. Physics of Fluids A, 1993, 5(12): 3186-3196.

[6] Gamet L, et al. Compact finite difference schemes on non-uniform meshes application to direct numerical simulations of compressible flows[J]. International Journal for Numerical Methods in Fluids, 1999(29): 159-191.

[7] Lenormand E, et al. Subgrid-scale models for large-eddy simulations of compressible wall bounded flows[J]. AIAA Journal, 2000, 38(8): 1340-1350.

[8] Mossi M, Sagaut P. Numerical investigation of fully developed channel flow using shock-capturing schemes[J]. Computers & Fluids, 2003(32): 249-274.

[9] 李新亮, 馬延文, 傅德薰. 可壓槽道湍流的直接數值模擬及標度律分析[J].中國科學(A輯), 2001, 31(2): 153-164.

Li Xinliang, Ma Yanwen, Fu Dexun. Direct numerical simulation of compressible channel turbulence and scaling law analysis[J]. Science in China (Series A), 2001, 31(2): 153-164.

[10] Coleman G N, Kim J, Moser R D. A numerical study of turbulent supersonic isothermal-wall channel flow[J]. Journal of Fluid Mechanics, 1995(305): 159-183.

[11] Lenormand E, Sagaut P, Ta Phuoc L. Large eddy simulation of subsonic and supersonic channel flow at moderate Reynolds number[J]. International Journal for Numerical Methods in Fluids, 2000(32): 369-406.

[12] 傅德薰, 馬延文. 高精度差分格式及多尺度流場特性的數值模擬[J]. 空氣動力學學報, 1998, 16(1): 24-35.

Fu Dexun, Ma Yanwen. High order accurate schemes and numerical simulation of multi-scale structures in complex flow fields[J]. Acta Aerodynamica Sinica, 1998, 16(1): 24-35.

[13] Lele S K. Compact finite difference schemes with spectral-like resolution[J]. Journal of Computational Physics, 1992(103): 16-42.

[14] Niederschulte M A, Adrian R J, Hanratty T J. Measurements of turbulent flow in a channel at low Reynolds numbers[J]. Experiments in Fluids, 1990, 9(4): 222-230.

[15] Kim J, Moin P, Moser R. Turbulence statistics in fully developed channel flow at low Reynolds number[J]. Journal of Fluid Mechanics, 1987(177): 133-166.

[16] Dean R B. Reynolds number dependence of skin friction and other bulk flow variables in two-dimensional rectangular duct flow[J]. Journal of Fluids Engineering, 1978, 100(2): 215-223.

Large Eddy Simulation of Compressible Channel Flow

Wang Suo-zhu, Yang Tian-peng, Yang Yi-feng, Qin Xu-guo, Su Wei

(Beijing Institute of Space Long March Vehicle, Beijing, 100076)

The compressible turbulent boundary layer is a basic scientific problem for the research of high-speed aircraft, and the understanding of its flow mechanism is of great significance for the aerodynamic and thermal design of high-speed aircraft. Large eddy simulation (LES) of compressible channel flow is performed. The LES method for compressible flow is established based on Favré filter, and a dynamic mixed subgrid-scale (SGS) model for compressible flow is formed by extending the mixed SGS model to the dynamic form. Through the analysis of SGS model, grid density and numerical scheme, the LES method established in the present work can obtain reliable results with fewer grids. The numerical results also show the transition process of the compressible channel flow, and capture the evolution process of the large-scale vortices and the near-wall velocity streaks.

compressible turbulence; channel flow; large eddy simulation

2097-1974(2023)01-0112-06

10.7654/j.issn.2097-1974.20230122

V211.3

A

2022-11-11；

2022-12-31

王鎖柱（1984-），男，博士，高級工程師，主要研究方向為再入飛行器空氣動力學設計。

楊天鵬（1993-），男，博士，工程師，主要研究方向為再入飛行器空氣動力學設計。

楊依峰（1990-），男，高級工程師，主要研究方向為再入飛行器空氣動力學設計。

秦緒國（1982-），男，博士，研究員，主要研究方向為再入飛行器總體設計。

蘇 偉（1979-），男，博士，研究員，主要研究方向為再入飛行器空氣動力學設計。