基于數值模擬的流場附面層邊緣識別方法

王 磊 張堃元 司江濤 劉凱禮 蔡北京 楊心宇

(1.上海飛機設計研究院，上海 201210；2.南京航空航天大學，南京 210016)

0 引言

附面層是壁面附近流動因黏性而受阻滯的一部分區域，其邊緣通常取在速度達到主流0.99倍的位置。對附面層的研究是進氣道、噴管、機翼等內外流流場設計中的重要內容，例如對無黏設計結果進行黏性修正、分析附面層對主流的影響、判斷流動是否有分離的趨勢等。獲取附面層厚度以及附面層內參數分布是進行這些研究的基礎。

工程設計中有一些簡單的經驗公式可用于估算附面層的發展。例如，文獻[1-4]計算噴管附面層位移厚度δ1采用發展長度x的線性函數：

δ1=xtanβ

(1)

其中，角度β與噴管出口馬赫數有關；文獻[5-6]根據公式(2)計算高超聲速進氣道附面層位移厚度：

δ1=Ax+Bxe-x

(2)

其中，A、B為根據經驗確定的系數；文獻[7]中應用了一種形式更復雜的經驗公式。這些公式計算量小、使用方便，但顯然其計算精度以及能夠獲取的信息是比較有限的。

求解附面層方程是較為傳統的附面層近似計算方法。其中多數是對附面層動量積分方程補充經驗關系式后求解，例如文獻[8-11]在進氣道、噴管、葉柵等流道設計的黏性修正中應用的方法。也有些直接針對附面層偏微分方程組進行數值計算，例如Cebeci[12]、杜國樑[13]所提出的方法，其中Cebeci提供的計算程序在許多進氣道、噴管設計中得到應用[14-18]。

隨著CFD技術的飛速發展，基于黏性N-S方程組的數值模擬已得到普遍應用，與上述兩種方法相比，其精度更高，適應性更廣。但N-S方程組求解中一般不特意區分附面層與主流，這使得如何在流場數據中準確認定附面層邊緣成了新的問題。按照通常的定義，附面層邊緣取在速度達到主流0.99倍的位置，而在較復雜的流場中主流本身并不均勻，因此需要根據經驗人為給定主流速度的數值[19-21]，有較大的主觀性并且難以實現自動化識別；有文獻中根據附面層和主流中法向速度梯度大小的差異判斷附面層邊緣[22]，也存在類似的問題。

針對在流場數值計算結果中識別附面層的困難之處，本文根據附面層的概念，提出以排除了附面層干擾的“參考主流流場”中的速度代替實際流動主流速度確定附面層邊緣的方法，以典型的超聲速壓縮型面為例對該識別方法進行了驗證，并根據識別結果考察了應用附面層修正的設計效果。

1 附面層邊緣識別方法

針對在流場數值計算結果中識別附面層的困難之處，本文根據附面層的概念，提出以排除了附面層干擾的“參考主流流場”中的速度代替實際流動主流速度確定附面層邊緣的方法，并以超聲速壓縮型面為例對該識別方法進行了驗證。

如引言所述，根據數值模擬結果識別附面層困難之處在于復雜流場中主流速度的確定，主流應是附面層外、沒有受到壁面與流體之間剪切力影響的流動。滿足這種要求的流場通過零剪切滑移壁面條件下的數值模擬可以獲得，這里稱之為“參考主流”，其中流動沒有附面層，完全排除了壁面與流體之間剪切力的影響。

在流場中黏性相互作用較弱時(對其強弱的判斷參見文獻[23]的說明)，可以認為實際流動中主流流動與這種參考主流相同。因此根據附面層的定義判斷其邊緣時，可以用這種參考主流中對應位置的速度代替所需的實際流場主流速度，這種方法避免了人為指定主流速度的主觀性，因而更為科學合理，也有利于識別過程的自動化。

圖1為該方法識別過程的示意圖。在壁面上的指定位置(例如圖1中x)建立沿壁面法向的單位向量n和沿壁面切向的單位向量s，沿向量n方向新建一系列網格點，依次驗證各點參數是否滿足公式(3)：

圖1 附面層邊緣識別過程示意圖

(V-0.99Vref)·s≥0

(3)

式中，Vref為上述參考主流中該網格點的速度向量，V為實際流動中的速度向量(可以是數值計算或實驗測量的結果)，第一個滿足該式的點即判定為該截面附面層邊緣。對于壁面所有位置重復該過程，能夠獲取一系列邊緣點，依次連接就得到了整個壁面上的附面層邊緣曲線。

值得說明的是，對于主流湍流度較小并且不存在強剪切的流動，也可以用無黏條件下的模擬結果近似代替滑移壁面條件下的計算作為參考主流，這有利于減小計算量。

本文在商業軟件ANSYS FLUENT的Scheme語言環境中編制了能夠自動完成整個計算和識別過程的程序。

2 流場數值模擬方法及算例校驗

本文研究中采用軟件ANSYS FLUENT進行數值模擬，首先利用文獻[24-25]中超聲速壓縮拐角的實驗結果考核了所采用的軟件和算法對超聲速流動中附面層計算的可靠性。

圖2為該實驗模型的示意圖。風洞實驗段尺寸203 mm×203 mm，來流馬赫數2.87，總壓0.68 MPa，總溫282 K，實驗模型為8°的斜楔，模型寬度152 mm，貼于風洞底面安裝，斜楔拐角距離噴管出口1.2 m。實驗采集了壓縮拐角前后若干截面上的參數分布，本文采用了拐角前25.4 mm(記為截面1)和拐角后66.04 mm(記為截面2)兩個垂直于壁面的截面上時均速度分布。實驗的其他數據見文獻[25]。

圖2 實驗模型示意圖[24]

針對該實驗進行了二維數值計算。計算模型參考了文獻[26-27]的做法，在截面1前補充長度為2.0 m的附面層發展平板以保證截面1附面層速度分布與實驗相同，圖3為計算模型示意圖。計算的總網格量約400 000，壁面y+約25。進口邊界條件為壓力遠場，出口為壓力出口，壁面為絕熱無滑移固壁邊界。

圖3 計算模型示意圖[24]

采用基于密度的耦合隱式算法求解RANS方程，對流通量采用基于二階迎風格式的AUSM方法計算。湍流模型選用標準k-ε模型，采用二階迎風格式離散，近壁采用標準壁面函數處理。比熱比取定值1.4，黏性系數采用Sutherland公式計算。以來流條件初始化開始計算，當各項殘差不再下降并且進出口流量相對誤差低于10-3時完成計算。

圖4顯示了截面1和截面2速度分布的實驗與計算結果，圖中η為沿壁面法向的坐標，uξ為平行于壁面方向的速度大小。其中截面1計算與實驗最大誤差為5.5%，均方根誤差0.35%，表明拐角前來流滿足實驗條件；截面2計算結果與實驗相吻合，速度最大誤差為4.4%，均方根誤差0.27%，表明所采用的計算方法是可靠的，后續研究中均采用類似的算法進行計算。

圖4 兩個截面上速度分布的實驗和計算結果

3 附面層識別方法的驗證

取來流條件為：M=6，p=2 511 Pa，T=221.6 K，在無黏條件下設計了以下兩種超聲速壓縮型面：

1)三級斜楔壓縮型面：三個壓縮角分別為5.07°、5.71°和6.47°，無黏情況下三道激波匯聚點距離型面前緣高度為100 mm，壁面長度為433 mm；

2)彎曲壓縮型面：根據壓力分布反設計得到，設計方法參見文獻[28]，無黏情況下激波高度100 mm，壁面長度L=394.5 mm，壁面壓力分布為分段函數：

(4)

其中,下標i為分段節點序號，p1/p=1.67，式中其他參數取值見表1。

表1 壁面壓力分布函數中的參數取值

對這兩個壓縮型面進行數值模擬得到流場參數，計算中來流湍流度為1%，來流湍流黏性與分子黏性之比為1。之后對兩個壓縮型面進行無黏計算作為參考主流，應用第1節所述識別算法提取了附面層邊緣。

圖5顯示了兩個壓縮流場的馬赫數等值線以及提取出的附面層邊緣曲線。從圖中可見，除有激波穿過的局部范圍外，所提取的邊緣曲線光滑連續，均位于速度急劇變化區域的外緣，與附面層的物理意義相符。

(a)三級斜楔壓縮

圖6為兩個算例出口截面近壁區域的速度和馬赫數分布曲線。其中三級斜楔壓縮出口截面主流參數均勻，按照0.99倍的定義，附面層厚度應為6.76 mm；而本算法識別的結果為7.04 mm(參見圖6(a))，兩者相對誤差僅4.1%，表明識別的結果是比較準確的。

從圖6(b)可見，彎曲壓縮型面出口截面主流不均勻，馬赫數和速度分布在遠離壁面時均有明顯增加。按照傳統的做法，需要人為選定主流速度的數值，不僅存在較大的主觀性，也難以實現自動化的識別。而本算法所識別的附面層邊緣位于速度梯度有明顯變化的位置，符合附面層的概念。

(a)三級斜楔壓縮

4 應用附面層修正的型面設計

在工程應用中設計進氣道、噴管等流道型面時，常采用先進行無黏設計、再修正附面層影響的方法。在進行修正設計時，通過本文所述的附面層邊緣識別方法，能夠更準確、方便地獲得各截面附面層邊緣位置及附面層內參數分布，之后容易根據定義計算其位移厚度δ1：

(5)

其中，δ為附面層厚度，ρ、uξ分別為密度和沿壁面切向的速度分量，下標δ表示附面層邊緣位置的參數。對于無黏設計的壁面坐標向外移動該距離可修正附面層對流量的影響。

應用此方法對上述彎曲壓縮型面進行了修正。圖7為修正前后壁面形狀、激波形狀的對比，圖8為修正前后壁面壓力分布及與設計值的對比。可見修正后激波高度與特征線法設計結果的相對誤差大約從2.0%降低到0.3%，壁面壓力與給定壓力分布之間的誤差也明顯降低，末端壓力的相對誤差大小|ε|從6.6%降低至2.3%。

圖7 修正前后的壁面形狀和激波形狀

圖8 修正前后壁面壓力分布及與設計值間的誤差

5 結論

1)研究了依據數值模擬結果識別附面層的方法，與依靠經驗公式、附面層方程的近似計算相比有更高的精度和更廣的適應性；

2)提出了用“參考主流”代替實際主流識別流場附面層邊緣的方法：以零剪切力滑移壁面邊界條件下數值模擬得到的不受附面層干擾的流場作為參考主流，在根據附面層定義確定其邊緣時，以該參考主流流場中的速度代替實際的主流速度進行判斷。這種方法避免了傳統做法中人為指定主流速度的主觀性，因而更為科學合理，也有利于識別過程的自動化；

3)對兩個超聲速壓縮流場附面層識別的測試表明，該識別方法獲得的附面層邊緣比較準確，其中斜楔壓縮出口截面上附面層厚度與采用實際主流速度判斷得到的厚度相對誤差僅4.1%；

4)根據本算法識別的附面層參數計算位移厚度對彎曲壓縮型面進行黏性修正后，激波高度與無黏設計結果之間的相對誤差從修正前的2.0%降低至0.3%，壁面末端壓力相對誤差從6.6%降低至2.3%，表明該方法可取得良好的效果。