可壓縮邊界層的入口合成湍流生成方法1)

王 天 孫 東 郭啟龍 李 辰 袁先旭 李 博

* (中國空氣動力研究與發展中心,空天飛行空氣動力科學與技術全國重點實驗室,四川綿陽 621000)

引言

湍流邊界層在航空航天工程中普遍存在,對飛行器的氣動特性有重要影響,因此實現對湍流邊界層的準確模擬具有迫切的實際需求.湍流數值模擬方法主要包括[1]直接數值模擬(direct numerical simulation,DNS)、大渦模擬(large eddy simulation,LES) 和雷諾平均Navier-Stokes 方程(Reynolds averaged Navier-Stokes,RANS)模擬.DNS 能夠解析流場中的全部尺度,常被用于機理分析,但其對計算資源要求較高,無法應用于高雷諾數復雜外形中的工程計算.LES 計算量較DNS 更小且能夠解析流場中的精細流動結構,但在近壁區域的網格需求仍十分巨大,無法滿足當前的工程應用.RANS 方法通過對雷諾平均方程未封閉項進行建模大大降低計算量,廣泛應用于當前的工程計算當中.但RANS 方法無法準確模擬大分離流動,也無法獲得精細的非定常流動特征.

過去20 年來,融合了LES 與RANS 優勢的混合RANS-LES 方法[2-4](hybrid RANS-LES method,HRLM)為發展滿足工程需求的湍流模擬手段提供了新思路[5].根據RANS 和LES 混合方式可分為分區混合方法和非分區混合方法.在分區混合方法中,需要人為確定交界面,并在兩側分別使用RANS 和LES.從RANS 向LES 的過渡過程中,由于缺少LES 所需的湍流脈動,流場中會產生一段適應區域使流場恢復至真實湍流.為縮短此區域范圍并節省計算資源,在交界面上添加合理的非定常脈動是十分必要的.

針對如何添加合理的湍流脈動這一問題,Tabor等[6]、Wu[7]和Dhamankar 等[8]針對這方面的研究進展做了全面詳盡的綜述.入口處的脈動生成方法被分為3 大類: 預先模擬類(precursor simulation)、回收調節類(recyling-rescaling method,RRM)、合成湍流法(synthetic turbulence,ST).3 類方法各有優缺點及其適用場合.從工程應用的角度來看,ST 類方法具有兩方面優勢: 一是其計算代價相對較小;二是其所需的輸入信息常規RANS 方法均可提供.因此ST 類方法得到了廣泛的關注.

ST 類方法本質上是從湍流的低階統計信息中重構高階統計信息,可以通過多種途徑實現.其中合成湍流生成器[9-11](synthetic turbulence generator,STG)、合成渦方法[12-13](synthetic eddy method,SEM) 和數字濾波法[14-15](digital filter method,DFM)是3 種應用較多的方式.這些方法的性能主要依賴于輸入的湍流低階統計量,如湍動能和長度尺度,其中長度尺度具有較大的不確定性,不同RANS 模型給出的近似分布可能存在明顯差別,Guo 等[11]討論了不同長度尺度分布對ST 性能的影響,并給出了一種基于一方程S-A 模型的長度尺度構造方式.另一方面,上述ST 方法添加的湍流脈動一般僅含速度脈動,在可壓縮邊界層中,除了包含速度脈動,還包含溫度、密度、壓力等熱力學變量的脈動.在高馬赫數湍流邊界層的HRLM 模擬中,也可以忽略熱力學量脈動,直接使用ST 方法在入口/交界面上添加速度脈動,但可能會降低合成邊界層脈動恢復到物理真實脈動的速度.對于如何在湍流入口中添加熱力學脈動這一問題,Martin[16]在開展可壓縮壁湍流直接數值模擬的過程中,基于強雷諾比擬(strong reynolds analogy,SRA)給出了流場初始的溫度脈動及密度脈動.SRA 是Morkovin[17]在研究可壓縮湍流中提出的直接反映速度脈動和熱力學量脈動之間關系的理論.此理論使得湍流入口的熱力學脈動量可直接由速度脈動得到.Adler 等[18]結合DFM 和SRA 給出了湍流入口的速度脈動和熱力學脈動.但是原始SRA 關系式被證明并不嚴格成立[19-20],Gaviglio[21]和Huang 等[22]提出了新的改進的SRA關系式.這些關系式均可用來構建入口的熱力學脈動.

本文主要開展兩方面工作: 一是在不可壓縮和可壓縮的典型壁湍流算例(槽道湍流和邊界層湍流)中,對比分析STG,SEM 和DFM 3 種合成湍流方法作為湍流入口對流場發展過程的影響;二是針對可壓縮壁湍流的模擬,探討了若干強雷諾比擬方法所得到的熱力學量脈動對流場發展的影響.

1 合成湍流入口方法及計算模型簡介

合成湍流方法通過局部的平均量和湍流統計特征量構建湍流脈動,操作簡單,被廣泛應用于湍流入口.

速度ui(r,t) 可以被分解為時間平均速度(r) 和脈動速度(r,t),如式(1)所示,r為位置矢量.合成湍流的目的是構建合理的脈動速度

在湍流邊界層中脈動速度具有各向異性,Davidson[23]提出了通過求解雷諾應力張量的特征值實現各向異性的方法,對雷諾應力張量R進行Cholesky 分解,即R=ATA

各向異性脈動速度通過如下方式獲得

1.1 合成湍流生成器

本文采用的STG 方法為文獻[11]中提出.首先通過一系列時空Fourier 模態加權疊加得到輔助矢量

其中上標“n”表示第n個模態的參數,q是歸一化幅值,由RANS 結果計算得到,kn=kn·dn為三維的波數矢量,kn表示波數矢量的模,dn為均勻分布于單位半徑的球面上的隨機方向矢量,σn為隨機生成的、與dn垂直的單位矢量,φn是均勻分布于[0,2π]區間的隨機相位.所有以上這些隨機量隨n變化,且在計算過程中不隨時間變化.具體細節及構造方法詳情可參見文獻[10].

1.2 數字濾波法

數字濾波法最早由Klein 等[24]提出,此方法利用濾波的數學性質,設計滿足在隨機波動中建立高斯兩點相關函數的數字濾波器,構建滿足所需長度尺度的隨機速度分布.之后Xie 等[25]采用二維濾波改進了原始的方法,大大節省了計算資源.Adler 等[18]提出了一種改進的數字濾波法能夠滿足復雜外形.本文采用文獻[18]給出的構造方式,其速度脈動采用式(3)構建,輔助矢量的構建方法如下

式中 ηl(j,k) 為定義在入口或交界面平面上的濾波隨機向量,(j,k) 表示平面上的網格計算坐標,It為積分時間尺度.

式中,下標0 為選定的流場位置,式(6)即以選定位置為中心進行濾波,α 為歸一化常數,Ij和Ik為各方向上的積分長度尺度,r為隨 (j,k) 均勻分布的隨機數序列.

1.3 合成渦方法

Jarrin 等[12]通過在入口平面引入人工渦的方法建立速度脈動場,每個渦由特定形狀函數表示,描述其空間和時間特征.本文采用文獻[13]中速度脈動構建方法

其中,VB為旋渦所占體積,σ 為模型的長度尺度,由RANS 結果得到.其他具體參數可參見文獻[13].

1.4 數值方法和模型

本文對可壓縮Navier-Stokes 方程進行求解,采用基于S-A 模型的IDDES 進行模擬,詳情見文獻[26].其中無黏項采用文獻[27]中的hyWENOUW4 格式,黏性項離散使用4 階中心差分格式,時間推進采用隱式雙時間步法.本文算例中均在入口添加人工合成湍流脈動.

2 合成湍流入口條件對比分析

本節分別采用STG,DFM 和SEM 方法生成湍流脈動應用于不可壓縮槽道湍流和可壓縮壁湍流的入口條件中,對比分析了其對流場壁面摩阻、流場結構、雷諾應力等隨流場發展的影響.

2.1 對不可壓縮槽道湍流發展影響

槽道湍流被廣泛用于對湍流入口條件的生成方法進行測試[10,28-29].本節分別使用STG,SEM 和DFM 方法在入口添加湍流速度脈動對充分發展槽道湍流進行模擬.來流馬赫數為0.2,基于槽道半寬h和壁面摩擦速度uτ的雷諾數Reτ=395,計算域范圍為x×y×z=32.5h×2h×3.4h.流向(x方向)、法向(y方向)和展向(z方向)網格數分別為Nx×Ny×Nz=469×109×139,其中第1 層網格高度ymin=0.02h,流向網格和展向網格均勻分布Δx=0.07h,Δz=0.025h,對應的≈0.79,Δx+≈27.65,Δz+≈9.875.槽道展向采用周期性邊界條件,壁面法向邊界采用等溫無滑移條件.

圖1 展示了入口截面(x/h=0)上3 種合成湍流方法下入口的展向脈動速度云圖.STG 和DFM 均添加了大量的脈動區域,區域范圍沿壁面法向逐漸變大;在接近壁面處,DFM 添加的脈動范圍呈現法向短展向長的細長形;而SEM 添加的脈動數量較少,且在相同法向距離下脈動尺度明顯大于STG 和DFM.

圖1 入口處展向速度云圖Fig.1 Spanwise velocity contour of inflow

通常采用時均摩阻沿流向的變化來反映入口合成湍流的品質,定義恢復距離為時均摩阻完全恢復至參考值5% 以內的位置距入口的流向長度.圖2展示了STG,DFM 和SEM 3 種合成湍流入口條件下時均摩阻分布.RANS 能夠準確模擬此構型的流場[30],因此圖中以RANS 模型計算結果作為參考值.對比結果表明,在添加合成湍流后,摩阻均在入口附近的一段區域內明顯低于參考值.其中,DFM 和SEM 兩種方法得到的恢復距離約為15h,STG 得到的恢復距離約為7h,優于其他兩種方法.

圖2 不同合成湍流入口下摩阻對比Fig.2 Comparison of frictions of different inlet boundary conditions

圖3 給出了不同合成湍流入口條件下半槽道流場的瞬時Q等值面(Q=0.05),用于顯示旋渦結構,等值面采用流向速度著色.結果顯示,STG 方法和DFM 方法添加湍流脈動后,入口后出現類似的大尺度流向相干結構,其中DFM 流場中此結構區域較大,導致恢復距離明顯長于STG.SEM 添加湍流脈動后無大尺度流向相干結構.圖4 給出了y/h=0.2位置x-z截面的渦量云圖.對比STG 和DFM 的結果,在流場入口后均存在大尺度結構,但STG 的恢復距離明顯更短.SEM 下流場結構隨流向變化并不明顯.

圖3 流向速度著色的Q 等值面(Q=0.05)Fig.3 Isosurfaces of the Q-criterion coloured by streamwise velocity(Q=0.05)

圖4 槽道y/h=0.2 橫截面渦量云圖Fig.4 Contours of vorticity,over the x-z plane at y/h=0.2 in channel flow

為定量描述流動的發展過程,圖5～圖7 給出了3 種湍流入口下流向x/h=0,3,6,12,20 位置處的基于摩擦速度的雷諾應力.其中選擇周期邊界條件的IDDES 結果及Poletto 等[31]計算結果代表完全發展湍流值作為參考.圖5 展示了R12隨流向的發展過程,可以發現,STG 方法下切應力向參考值恢復得更快.其中在x/h=6 處,STG 下已基本恢復,而DFM 和SEM 下的值仍與參考值相差較大.在x/h=12 后,各方法下R12已基本恢復至參考值,但繼續隨流向發展DFM 下R12與參考值仍有一定差異.對比圖6 中R11發展過程,STG 及SEM 下均能迅速恢復至參考值,而DFM 下恢復速度較慢,在x/h=12 后才逐漸與參考值靠近.對于圖7 中R22發展變化曲線,3 種入口條件下恢復均較慢,在x/h=12處數值均與參考值有一定差距.

圖5 不同合成湍流入口下 R12 隨流向變化Fig.5 R12 development along streamwise using different inlet conditions

圖6 不同合成湍流入口下 R11 隨流向變化Fig.6 R11 development along streamwise using different inlet conditions

圖7 不同合成湍流入口下 R22 隨流向變化Fig.7 R22 development along streamwise using different inlet conditions

2.2 對可壓縮壁湍流發展影響

針對可壓縮算例,本節選擇Ma=2.5 零壓力梯度平板算例開展研究,流向(x方向)、法向(y方向) 和展向(z方向) 網格數為Nx×Ny×Nz=374×121×59,其中法向第1 層網格≈0.38,流向網格和展向網格均勻分布 Δx+≈24,Δz+≈13,出口處網格稀疏處理,使用海綿層[32]來消除出口邊界的虛假反射波.表1 中 θ 為入口動量厚度.

表1 來流條件Table 1 Inlet conditions

圖8 展示了超聲速平板入口的展向速度脈動,δ0為入口邊界層厚度.可以直觀地觀察不同合成湍流方法在入口處添加的脈動.可以發現,STG 和DFM 流場中接近壁面處均存在多個小范圍速度脈動區,DFM 中湍流脈動區域更多呈現展向細長狀,且尺度較小集中于壁面附近.SEM 方法添加湍流脈動后,入口處存在數個大范圍速度脈動.

圖9 展示了STG,DFM 和SEM 3 種方法下壁面摩阻沿流向的發展曲線,其中選擇Van Driest[33]的經驗公式作為對比.圖中,STG 和SEM 下摩阻發展趨勢基本一致,恢復距離約為 10δ0,而DFM 摩阻一直偏低,且恢復緩慢,發展至x/δ0=20 處仍有一些差距.總的來說,Ma=2.5 超聲速平板流中STG 和SEM 入口條件下恢復距離相比較短,DFM 恢復距離最長且精度最差.

圖9 不同合成湍流入口摩阻對比Fig.9 Comparison of frictions of different inlet boundary conditions

圖10 展示了邊界層內流場的Q等值圖(Q=0.2),其中采用流向速度著色.圖中觀察到STG 與DFM流場發展過程相似,入口處存在流向相干結構,但STG中此結構明顯短于DFM 中.SEM 下流場逐漸由小尺度結構發展為較大尺度的湍流結構.圖11 給出了y/δ0=0.2橫截面的渦量云圖,與上面描述基本一致.

圖10 邊界層內Q 值等值面(Q=0.2)Fig.10 Isosurfaces of the Q-criterion in boundary layer (Q=0.2)

圖11 湍流邊界層 y/δ0=0.2 橫截面渦量云圖Fig.11 Contours of vorticity,over the x-z plane at y/δ0=0.2 in turbulent boundary layer

本節分別選擇流向x/δ0=0,3,6,12,20 位置來對比不同合成湍流方法下基于黏性尺度的雷諾應力隨流向的發展變化過程,如圖12～圖14 所示.其中選擇Pirozzoli 等[34]的DNS 結果作為參考.對于而言,SEM 入口條件下迅速恢復至參考值附近,STG 中流場不同流向位置處曲線變化并不大,但與參考值有一定差距,而DFM曲線隨流向變化劇烈且誤差更大;隨流向變化中可以發現,STG 和SEM 下曲線可迅速恢復至參考值,且SEM 精度更高;而隨流向變化圖中STG方法下精度最高.

圖12 不同合成湍流入口下 隨流向變化Fig.12 development along streamwise using different inlet conditions

圖13 不同合成湍流入口下 隨流向變化Fig.13 development along streamwise using different inlet conditions

3 添加熱力學脈動入口條件對比分析

本節選擇Ma=2.5,5.0 工況下開展入口添加熱力學脈動對流場發展影響研究.其中Ma=2.5 計算條件與2.2 節中相同.在Ma=5.0 算例中,流向(x方向)、法向(y方向) 和展向(z方向) 網格數為Nx×Ny×Nz=336×181×69,其中法向第1 層網格≈0.46,流向網格和展向網格均勻分布 Δx+≈36.8,Δz+≈18.4.來流條件如表2 所示.

表2 來流條件Table 2 Inlet conditions

3.1 熱力學脈動添加方法

使用式(3)僅能給出速度脈動,但對于高馬赫數邊界層,熱力學量的脈動也是不可忽略的一部分.通常的做法是使用各種強雷諾比擬(strong Reynolds analogy,SRA)來給定溫度脈動的均方根值.在絕熱壁面條件下,溫度脈動的均方根可以表示為

式中T表示溫度,“～”表示Favre 平均量.但是對于等溫壁面,式(11)給出的統計結果并不完全正確,又有學者提出了修正的雷諾比擬[35],其統一形式可以寫為

本文所采用的原始強雷諾比擬以及Gaviglio[21]和Huang 等[22]修正后的參數a和c的取值參見表3,分別稱為SRA,GSRA 和HSRA.其中HSRA中c=Prt,如下式所示.但在入口處Prt需要取定值,Zhang 等[35]提到Prt參考值為 0.7 ～0.9.本文為使GSRA 及HSRA 添加熱力學脈動入口條件后流場發展過程對比明顯,在入口中取Prt為0.7.

表3 入口條件中強雷諾比擬系數表Table 3 SRA’s coefficients of inlet conditions

在由式(12)得到溫度脈動后,通過完全氣體狀態公式和零壓力脈動假設可以得到密度脈動,如下式所示

3.2 入口添加熱力學脈動對流場發展影響

綜合第2 節中的對比結果,STG 相對其他方法表現出一定的優勢.因此為對比可壓縮條件下不同熱力學脈動生成方法對流場發展的影響,本文采用STG 方法添加速度脈動,并通過不同強雷諾比擬在入口處添加熱力學脈動.本節分別模擬了Ma=2.5和Ma=5.0 條件下的湍流邊界層,后續采用“case 0”代表入口處不添加熱力學脈動,分別采用“case 1”,“case 2”和“case3”代表采用SRA,GSRA 和HSRA 添加熱力學脈動.圖15 給出了不同馬赫數下得到的摩阻分布.從摩阻分布中發現,不同熱力學脈動入口條件下恢復距離基本一致.圖16 給出了完全湍流區內基于黏性尺度的雷諾應力數值不同入口條件下結果差距較小,且均能取得合理的結果.

圖15 不同入口條件下摩阻曲線Fig.15 Comparison of frictions of different inlet boundary conditions

圖16 邊界層模擬得到的完全湍流區后的雷諾應力Fig.16 Reynolds stresses in fully-developed turbulence

下面將分析SRA,GSRA 和HSRA 添加熱力學脈動后的入口條件對流場發展的影響.式(13)中湍流Prandtl 數表征平均運動中動量交換和熱交換之比,在強雷諾比擬中此值近似為1.圖17 展示了Ma=2.5,5.0 下完全湍流區內不同入口條件下湍流Prandtl數沿法向的分布,其中線點圖表示Ma=2.5 條件,DNS 結果取自文獻[36].可以發現,在完全湍流區內,采用不同強雷諾比擬添加熱力學脈動作為入口條件后,湍流Prandtl 數變化趨勢基本一致,其數值沿法向從1 逐漸下降至0.7 附近,與Pirozzoli 等[37]和Duan 等[38]的結果也基本一致,證明了本文對于流場熱力學量模擬的準確性.

圖17 完全湍流區內湍流Prandtl 數法向分布Fig.17 Turbulent Prandtl number in fully-developed turbulence

在不同流動條件下(槽道和邊界層),不同壁溫條件下有大量DNS 數據[37-41]表明HSRA 相比于SRA 和GSRA 表現更好.Guarini 等[41]深入分析了HSRA 關系,指出HSRA 有效意味著湍流場中法向方向的動量和熱量輸運具有相似性.因此本文定義由HSRA 定義的公式來衡量流場熱力學量的發展過程

圖18 展示了兩種馬赫數下完全湍流區內f的數值,其中線點圖代表Ma=2.5.不同入口條件下不同馬赫數下曲線基本一致.在完全湍流區處f值在邊界層內沿法向逐漸降低至1 附近.

圖18 完全湍流區內 f 分布Fig.18 f distribution in fully-developed turbulence

圖19 給出了Ma=2.5,5.0 下x/δ0=2,6,10,14位置處的f的數值,其中線點圖代表Ma=2.5.隨流向發展,不同入口條件下f逐漸靠攏至1 附近.可以對比觀察到圖中case 0 曲線向完全湍流數值恢復得最慢,即在流場入口處不添加熱力學脈動后流場中熱力學量向完全湍流狀態恢復得更慢.在不同位置處,可以明顯觀察到case 2 和case 3 入口條件下流場熱力學量恢復得更快.在x/δ0=2,4 兩個位置處可以發現case 2 中恢復速度最快.因此采用GSRA,HSRA 在入口處添加熱力學脈動均能加快流場中熱力學量的恢復速度,且GSRA 效果稍好于HSRA.

圖19 Ma=2.5,5.0 下不同流向位置fFig.19 f distribution in various streamwise position of Ma=2.5,5.0

3 結論

本文分別以STG,DFM 和SEM 3 種方法添加入口湍流脈動開展了低馬赫數槽道湍流和超聲速湍流邊界層的數值模擬研究,通過對比分析流場中的湍流脈動、壁面摩阻、流場結構、雷諾應力和熱力學量等信息后,得到以下結論.

(1) 在不可壓縮和可壓縮(不添加熱力學量脈動)壁湍流問題模擬中,STG,DFM 和SEM 方法均能獲得充分發展的湍流結構.對比發現,在不可壓槽道湍流中,STG 能夠最快獲得合理的摩阻分布,流場結構與雷諾應力發展也有一定的優勢;在可壓縮湍流中,STG 與SEM 表現相近,均優于DFM 方法.

(2) 在馬赫數2.5 和5 下分別在入口處通過SRA,GSRA 和HSRA 添加熱力學脈動開展湍流邊界層數值模擬研究.結果顯示,是否添加熱力學脈動對于摩阻和雷諾應力發展影響較小,但對流場中的熱力學量有很大的影響.其中通過GSRA 添加熱力學脈動后流場熱力學量恢復最快.