李斯特 董 明

1 天津大學機械工程學院, 天津 300072

2 中國科學院力學研究所非線性力學國家重點實驗室, 北京 100190

1 引 言

對邊界層轉捩的有效預測是高超聲速飛行器設計中的基礎難題. 在環(huán)境擾動強度較低的“自然”狀態(tài)下, 存在三個與轉捩預測相關的基本問題(Kachanov 1994, 周恒和張涵信 2017). (1)感受性問題: 描述自由流中擾動如何進入邊界層并激發(fā)失穩(wěn)模態(tài)的過程; (2)線性失穩(wěn)問題: 描述失穩(wěn)模態(tài)是如何在邊界層中持續(xù)增長的; (3)轉捩判據(jù): 確定線性擾動累積到多大幅值時非線性轉捩發(fā)生. 工程中最常用的轉捩預測方法是eN法(Smith 1956, Cebeci et al. 1980, Su & Zhou 2009), 它主要基于線性穩(wěn)定性理論定量地描述機制(2), 而把非線性轉捩判據(jù)與由感受性機制決定的線性擾動初始幅值的比值當作人為參數(shù)處理. 該理論的優(yōu)點是簡單易用, 但也存在明顯不足: 首先, 該方法沒有考慮機制(1)與(3)的影響, 故人為參數(shù)的確定較為隨意; 更重要的是, 工程實際中飛行器表面很難做到完全光滑, 那里往往存在如粗糙元、縫隙、臺階等局部突變, 而它們對轉捩的影響是傳統(tǒng)的eN法無法考慮的. Schneider (2008a, 2008b)和Wheaton 等 (2010)的實驗表明, 較大尺度的粗糙元可以促進超聲速邊界層轉捩; 而Fujii (2006)和Fong 等 (2014)分別通過實驗與直接數(shù)值模擬(direct numerical simulation, DNS)發(fā)現(xiàn)二維粗糙元在一定條件下可以使轉捩推遲.Tang 等 (2015)通過風洞實驗方法研究了高超聲速邊界層中第二模態(tài)擾動在二維粗糙元附近的快速畸變過程. 為了揭示粗糙元對轉捩的影響機理, 需要進行更精細的理論研究工作. Wu 和Dong (2016)指出, 當失穩(wěn)模態(tài)的波長與局部突變的流向尺度相當時, 線性穩(wěn)定性理論與eN法均失效, 為此, 他們構建了一套新的理論框架?局部散射理論. 在該文章中, 他們關注的是亞聲速機制, 之后該理論被推廣到跨聲速機制(Dong 2020)以及超聲速與高超聲速機制(Zhao et al.2019, Liu et al. 2020, Dong et al. 2020, Dong & Zhao 2021). 實際上, 局部突變對轉捩的影響主要存在兩個機理, 分別是局部感受性(局部突變與自由流中的擾動相作用而激發(fā)額外的失穩(wěn)模態(tài))與線性模態(tài)的局部散射(局部突變與來流失穩(wěn)模態(tài)相作用而改變后者的幅值). 董明(2020)曾對局部散射理論在以上兩種機制下的研究進展進行了系統(tǒng)的綜述. 針對二維凹槽型局部突變, Dong和 Li (2021)還對這兩種機制下線性擾動的演化開展了細致的直接數(shù)值模擬, 并把結果與理論進行了系統(tǒng)地對比. 但是該模擬只局限在層流階段, 并未包含擾動的非線性轉捩. 本文擬設計一個直接數(shù)值模擬方案: 假設高超聲速邊界層中已有失穩(wěn)模態(tài)被激發(fā), 計算它們在一系列粗糙元的影響下觸發(fā)轉捩的過程, 并把數(shù)值結果與理論預測進行細致地對比, 以驗證局部散射理論在轉捩預測中的精度.

2 物理問題與數(shù)值方法

本文選取的物理模型為帶粗糙元的高超聲速平板邊界層, 如圖1所示. 來流條件與Fong 等(2014)、Zhao 等 (2019)和Dong 和 Zhao (2021)相同, 馬赫數(shù) Ma、 來流溫度與單位雷諾數(shù)分別為 5.92, 4 8.69 K和 1.32×107m?1. 假設壁面存在三個流向孤立、展向均勻分布、形狀相同的粗糙元, 它們到平板前緣的距離分別為 125.3 mm, 185 mm和 244.7 mm.

以第二個粗糙元中心位置原點建立笛卡爾坐標系. 通過可壓縮光滑平板邊界層的Blasius相似性解可以估算出該處的排移厚度本文選取這個值作為無量綱化的參考長度. 這樣, 三個方向的無量綱坐標為其中上標 * 為有量綱值. 三個粗糙元中心的無 量 綱 坐 標 分 別 為x1=?30,x2=0,x3=30. 令 每 個 粗 糙 元 的 形 狀 分 布 為h和d分別標識每個粗糙元的高度與寬度. 在本文的計算中,h=0.24,d=2.37 . 在計算中布置貼體網(wǎng)格, 如圖2所示, 并建立貼體坐標系 (ξ,η,ζ). 計算域取?42.65≤ξ≤230.65, 0≤η≤50, 0 ≤ζ≤23.64 , 計算網(wǎng)格數(shù)為 2129×151×128. 各物理量均以來流參量進行無量綱化, 進而定義雷諾數(shù)其中u∞和ν∞分別為來流流向速度與運動黏性系數(shù).

圖2

物理模型示意圖

控制方程為貼體坐標下的三維可壓縮Navier?Stokes方程組, 采用有限差分法對其離散. 在計算基本流時, 對流項采用五階WENO (weighted essentially non?oscillatory)差分格式, 黏性項采用四階中心差分格式, 時間推進采用LU?SGS (lower?upper symmetric Gauss?Seidel)方法; 在計算擾動時, 對流項采用五階迎風差分格式, 黏性項采用六階中心差分格式, 時間推進采用三階Runge?Kutta法, 具體離散過程見趙磊 (2017). 首先計算二維定?；玖? 即在圖1的三維計算域中選擇一個z=0的切片. 計算域入口給定為可壓縮Blasius解, 上邊界和出口采用外推邊界條件, 壁面為無滑移、無穿透、等溫的邊界條件, 其中無量綱壁面溫度Tw為 6.88. 二維基本流計算定常后, 把它展向均勻地布置在三維計算域中, 計算三維轉捩過程. 在入口引入 5個特征模態(tài)擾動,形式為其中下標i表 示第i個 擾動;分別表示密度、流向速度、法向速度、展向速度和溫度的特征函數(shù);ω為頻率,β為 展向波數(shù),A為擾動初始幅值,c.c.表示復共軛. 在三維非定常數(shù)值模擬中, 上下邊界條件與計算基本流時相同, 但出口改用嵌邊邊界條件, 以保證擾動無反射地傳出計算域. 展向采用周期邊界條件. 在給定的ω和β下 ,q?由可壓縮Orr-Sommerfield (O-S)方程給出. 在本文計算工況中, 選取的5個模態(tài)擾動的參數(shù)如表1所示.圖3展示了計算域入口處, 不同頻率擾動的增長率?αi, 圖中的“°”標識本文引入的擾動. 為定量刻畫粗糙元對轉捩的影響, 本文計算了光滑平板與帶粗糙壁面的兩種工況, 分別標記為Case 1與Case 2.

圖3

表 1 模態(tài)擾動參數(shù)

3 散射理論的應用檢驗

圖4為Case 2的平均流壓力的等值線云圖. 每個粗糙元對平均流的修正在無黏勢流區(qū)表現(xiàn)為壓縮?膨脹波系, 其馬赫角為 tan?1(Ma2?1)?1/2≈10°.圖5進一步展示了x=0附近的壓力等值線云圖和流線, 可以清楚地看到在粗糙元上游的高壓區(qū)與下游的低壓區(qū), 流線也在凸起附近明顯彎曲, 但在粗糙元前后均未出現(xiàn)分離, 其他粗糙元附近的流場也類似. 從定量上看, 該結果與Zhao 等 (2019)文中的圖4(c) 吻合很好. 在距離突變足夠遠的上下游, 平均流也趨近于未受擾動的Blasius解.

圖4

圖5

對于亞聲速邊界層中的黏性Tollmien?Schlichting (T-S)模態(tài), 粗糙元往往促進擾動的增長,從而使轉捩提前(Wu & Dong 2016); 但在高超聲速邊界層中, Fujii (2006)通過實驗發(fā)現(xiàn)了相反的結論. 加州大學洛杉磯分校的Zhong課題組的一系列數(shù)值工作(Duan et al. 2013; Fong et al. 2014,2015)表明, 粗糙元只會促進頻率低于同步頻率的擾動增長, 而對高頻擾動有抑制作用. Zhao等 (2019)采用求解諧波型線性化Navier?Stokes方程(harmonic linearized Navier?Stokes equation, HLNS) 在馬赫數(shù)為5.92工況下, 對透射系數(shù)與來流擾動的頻率、展向波數(shù)以及粗糙元大小和位置的關系進行了系統(tǒng)的數(shù)值計算. Dong 和Zhao (2021)采用高雷諾數(shù)漸近理論, 對該現(xiàn)象給出了深刻解釋. 本文的DNS可進一步展示該機制對非線性轉捩階段的影響.

對由DNS得到的流場進行傅里葉變換

其中qˇn,m(x,y)為 頻率為mω0、 展向波數(shù)為nβ0的 Fourier分量,ω0=0.2為 基頻,β0=0.7972為基波數(shù),T=2π/ω0為 時間周期,ZK=6π/β0為 展向計算域長度,q?(x,y,z,t)為瞬時擾動量.圖6給出了光滑壁情況下典型第二模態(tài)擾動(表1中的擾動3, 4和5, 它們所對應的 (m,n) 分別為(8, 0), (9, 0)和(11, 0))的傅里葉分量的幅值(傅里葉分量的模的法向最大值)與線性穩(wěn)定性理論預測結果的對比. 可以看出在x= 0以前, 幅值最大的擾動(11, 0)幾乎線性增長, 在幅值達到0.01的量級后進入非線性飽和階段. 而擾動(8, 0)和(9, 0)在下游可以累積到更大的幅值, 最大幅值的擾動(8, 0)的幅值峰值達到了0.042.圖7中的虛線展示了更多頻率下的傅里葉分量在更寬的流向區(qū)域內的演化規(guī)律. 高頻第二模態(tài)(m= 8, 9, 11)在非線性飽和后迅速衰減. 低頻的第一模態(tài)(m= 1, 2, 3)在前期增長率較低, 但它們可以持續(xù)增長到x=130左 右. 在x=180以后, 各階擾動的脈動強度趨于穩(wěn)定, 它們的幅值在0.01到0.02之間波動, 這說明轉捩即將完成. 這一過程與Dong 和 Luo(2007)對高超聲速尖錐邊界層自然轉捩的規(guī)律相同: 由于第二模態(tài)的線性增長率較高, 它們是轉捩前的主導擾動; 但是, 只有三維的低頻模態(tài)的幅值累積到非線性狀態(tài)的時候, 轉捩才能被觸發(fā);第二模態(tài)對第一模態(tài)的增長起到“催化”作用.

圖6

圖7

圖7中的實線表示Case 2的結果, 它們與Case 1的主要區(qū)別是: 在粗糙元附近, 擾動的幅值存在明顯的“跳躍”, 這進一步影響擾動在下游的演化幅值. 引入透射系數(shù)T(Wu & Dong 2016)來量化這一跳躍, 其定義為下游與上游漸近幅值之比, 這里的漸近幅值可通過Case 1的結果得到.對于每個頻率、展向波數(shù)組合 (m,n)的 擾動, 在三個粗糙元處的透射系數(shù)記作T1(m,n),T2(m,n)和T3(m,n).

Zhao 等 (2019)通過HLNS, 得到相同馬赫數(shù)和Reδ0?= 26 307工況下的透射系數(shù)隨頻率的分布. Dong 和 Zhao (2021)進一步用漸近理論討論了這一結果受Re的影響. 實際上, 本文中三個粗糙元所處位置對應的雷諾數(shù)Reδ0?分別為21 648, 26 307和30 228; 而從Dong 和 Zhao (2021)的結果來看, 這三個雷諾數(shù)所對應的T~ω曲線非常接近.圖8給出由DNS數(shù)據(jù)讀出的n=0, 不同頻率下的T1,T2與T3; 這里對頻率重新歸一化, 此處所采用的特征長度為粗糙元中心處對應的光滑平板邊界層排移厚度. 圖中還畫出HLNS的結果(選自Zhao et al. 2019)做對比, 它們總體吻合很好. 另外, 在HLNS計算中, 擾動被假設為線性的, 但在實際DNS中, 有些擾動已經達到了非線性狀態(tài)(幅值O(0.01)). 但是兩者仍然能給出較一致的透射系數(shù), 這說明非線性效應在局部散射過程中作用較弱. 而即使擾動發(fā)展到弱非線性階段時, 粗糙元對擾動的局部散射效應仍可近似被透射系數(shù)描述.

根據(jù)局部散射理論(Wu & Dong 2016, Zhao et al. 2019, Dong & Zhao 2021), 對于有k個局部突變的平板邊界層, 透射系數(shù)與轉捩閾值的變化 ΔN存在對應關系通過觀察DNS計算出的擾動幅值演化(圖7)發(fā)現(xiàn)ω=1.6,β=0 (m= 8,n= 0)的模態(tài)在轉捩前占主導地位, 這與基于線性穩(wěn)定性理論的eN方法給出的結論相同. 通過該模態(tài)在三個粗糙元處的透射系數(shù)可以算出轉捩因子的變化值 ΔN=0.3736. 在傳統(tǒng)的轉捩預測eN法中, 當最危險的擾動根據(jù)線性理論放大eN倍時, 認為轉捩發(fā)生. 而對于存在粗糙元的工況, 需要用由透射系數(shù)確定的 ΔN修正轉捩因子, 并認為當最危險的線性擾動首次放大 eN?ΔN時轉捩發(fā)生. 根據(jù)eN法, 按擾動由線性過渡到飽和狀態(tài)的幅值選取轉捩因子N. 對于本文工況, 選取的N值為3.23, 進而讀出Case 2曲線相對與Case 1曲線的前移量 Δx≈9, 這就是局部散射理論預測的轉捩前移量.

圖9給出了壁面摩阻系數(shù)沿流向的變化. 對于光滑壁面的工況(Case 1),Cf曲線先緩慢減小,再從x=113.5 處 開始劇烈升高. 對于Case 2,Cf曲線在粗糙元附近有較明顯下降, 但在不遠的下游向未受擾動的狀態(tài)恢復. 在本文模擬的工況中, 粗糙元促進了失穩(wěn)模態(tài)的增長, 這使得Cf曲線劇烈抬升的位置有所提前. 若把Cf曲線最低點看作轉捩的起始點, 則粗糙元的存在使轉捩提前了Δx≈7. 這與理論預測的轉捩前移量非常接近, 從而說明了局部散射理論的可靠性.

圖9

4 結論

為了驗證局部散射理論在高超聲速邊界層中轉捩預測的精度, 本文設計了一套直接數(shù)值模擬工況: 分別在光滑壁面與包含若干粗糙元的高超聲速邊界層中引入相同的初始失穩(wěn)模態(tài), 以定量考察粗糙元對轉捩的影響. 得到如下結論:

(1) 通過對比兩個工況下擾動傅里葉分量的演化規(guī)律, 定量刻畫了粗糙元對擾動演化的影響,該結果與局部散射理論的預測基本吻合. 同時, 數(shù)值結果進一步表明, 雖然局部散射理論假設擾動是線性的, 但其結果也可以近似應用于描述擾動在弱非線性階段的演化.

(2) 通過透射系數(shù)可以定量計算主導模態(tài)擾動的轉捩因子的變化值, 進而確定其對擾動幅值累積的影響. 由該機制得到的轉捩位置的變化量得到了直接數(shù)值模擬的驗證.

值得說明的是, 本文模擬的工況中, 由粗糙元引起的轉捩的提前量較小. 這是由于本文設計的粗糙元較少, 且每個粗糙元的高度較低. 在真實應用中, 更多更大尺寸的粗糙元可使轉捩前移更加明顯. 同時, 巧妙地設計粗糙元, 使得主導擾動的頻率盡量落在圖8的高頻區(qū), 是有可能使轉捩推遲的, 這將是作者下一步的工作.

附錄A

對本文直接數(shù)值模擬精度的驗證

由于風洞實驗中的背景擾動不易確定, 很難把數(shù)值結果與實驗結果做細致對比. 為了驗證本文計算的可靠性, 首先對計算網(wǎng)格進行驗證. 保留本文計算的網(wǎng)格寬度, 把流向計算域縮短為x∈[?42.6,107.9],計算網(wǎng)格數(shù)變?yōu)?201×151×128. 為了使流動在更短的距離轉捩,初始擾動((ω,β)=(0.5,0.7972), (1.0,0)和 (2.5,0)) 的幅值被增大為 0.01. 同時計算了網(wǎng)格數(shù)為 2401×151×128和1201×151×256 的流場, 并把三套網(wǎng)格下計算出的Cf曲線進行對比, 如圖A-1所示. 可以看出, 流向網(wǎng)格的加密對計算結果沒有影響, 而展向網(wǎng)格的加密使轉捩后期及湍流區(qū)的Cf值略大. 由于本文更關注的是轉捩前的局部散射作用與轉捩起始位置, 所以可以認為本文計算網(wǎng)格達到精度要求. 此外,選取Zhao 等 (2019)中圖9 (d)的工況, 重新計算了擾動的流向演化, 其擾動溫度如圖A-2所示. 兩條曲線精確吻合. 以上驗證說明了本文的計算結果是可靠的.

附圖 A-1

附圖 A-2

致 謝本文受到國家自然科學基金的資助 (U20B2003, 11772224).