壁面函數在超聲速湍流模擬中的應用

王新光,毛枚良,何琨,*,陳琦,萬釗

1. 中國空氣動力研究與發展中心 計算空氣動力研究所,綿陽 621000 2. 中國空氣動力研究與發展中心 空氣動力學國家重點實驗室,綿陽 621010

湍流數值模擬是當前計算流體力學(CFD)的核心問題之一,在現有高速大存儲計算機資源的支持下,涌現了大批關于直接數值模擬和大渦模擬的相關湍流研究成果,但對于復雜工程外形的湍流模擬,目前仍只能采用計算資源消耗較小的RANS方法。眾所周知,對于工程中常用的SST-模型,如需精確模擬壁面摩阻和熱流,通常要求在壁面網格密度滿足≈1,所需計算機資源往往難以滿足工程應用對效率的要求。為此,主流商業軟件中通常采用壁面函數,耦合不同的湍流模型來進行湍流模擬。

壁面函數從理論研究角度來看,是邊界層方程的近似解,其最初形式為對數律,刻畫了零壓力梯度不可壓附體邊界層流動的速度分布,之后得到了考慮壓力梯度和零剪應力情況下邊界層方程的漸近解,文獻[2-4]都相繼對標準壁面函數進行了修正研究。其中Nichols發展了流動可壓縮性和熱傳導的壁面函數,并將其應用到可壓縮湍流平板的數值求解中,賀旭照等將其與-兩方程模型耦合應用到可壓縮拐角的數值求解中。另外國內吳曉軍和肖紅林等分別在RANS模擬和大渦模擬簡單二維流動中采用壁面函數,Gao和Tao等分別通過數值實驗和風洞試驗對壁面函數的系數進行了壁面函數修正研究。但壁面函數在壓力梯度不可忽略的區域內是否仍然有效尚無定論,對于存在分離和再附的復雜流動依然存在困難。

文獻[15]通過DNS總結壁面附近湍流黏性系數分布,解析求解湍流邊界層方程,從而獲得壁面網格內的解析速度和溫度分布,稱為解析壁面函數,并將其推廣到可壓縮流動,通過簡單流動驗證了解析壁面函數對于在分離流動中具有良好的表現。

另一方面壁面函數僅作用于壁面網格,與全局流動的湍流模型存在不相容的問題,容易引起數值模擬結果對壁面附近第一層網格節點位置的強烈依賴性,使得數值預測的分離和再附區域精度變差。Tidriri和Knopp等將計算區域劃分為近壁區域和整體區域,在近壁區域通過子迭代的方式求解壁面剪切應力,在固壁面上始終滿足原始邊界條件,在近壁區域的外緣和整體RANS數值求解匹配一致。Kalitzin等通過不可壓縮平板邊界層數值計算,建立了摩擦速度和壁面第一個點的雷諾數之間的關系數據表,求解摩擦來得到壁面剪切應力,從而數值實現壁面函數在其不可壓縮笛卡爾RANS求解器中的應用。

目前大多數壁面函數研究主要基于平板邊界層,而工程實際的飛行器表面,通常是三維曲面,這使得一些復雜的壁面函數在實際應用中容易出現數值剛性的問題,該方法對于工程復雜外形的推廣應用仍有一定距離,這也是在一些復雜流動的模擬中又重新回歸到一些簡單的壁面函數的原因。

本文基于國家數值風洞氣動力計算平臺, 針對工程湍流模型SST-,開展適用的壁面函數研究,給出了具體的程序實現過程和邊界處理形式,通過典型二維算例完成了壁面函數模塊考核,最后基于內外流一體化的超聲速飛行器湍流流動的模擬,完成了模塊對復雜外流的工程應用考核。

1 壁面函數

對于充分發展的湍流邊界層,基于所謂的雙層模型,分為黏性底層和對數率層,如圖1所示。

黏性底層是一個緊貼固體壁面的極薄層,分子黏性切應力是主要特征因素,湍流切應力極小且可以忽略,所以流動似乎處于層流狀態,沿壁面法向呈線性分布,其區域為<。對數率層是一個距壁面較遠的流體層,流動處于完全湍流狀態,湍流渦黏性遠大于分子黏性,速度分布接近于對數率,位于從>到0.2倍邊界層厚度范圍內,在局部平衡的湍流邊界層內,可通過理論推導得到標準壁面函數,其具體形式為

(1)

(2)

式中:

=()12=

(3)

=()12=

(4)

圖1 湍流邊界層示意圖Fig.1 Sketch map of turbulence boundary layer

2 壁面函數在三維邊界層中的推廣應用

從平面推廣到三維曲面,就是把邊界層方程及其相關理論表達式轉換到貼體坐標系中,如圖2 所示。以當地離開壁面第1個計算單元處的切向速度作為流動方向,將三維流動進行局部二維化處理,在此基礎上采用湍流邊界層的壁面函數,得到壁面剪切應力。

圖2 三維壁面網格示意圖Fig.2 Sketch map of three-dimentional near-wall mesh

為了盡可能減少對源程序的修改,通過引入壁面虛擬湍流渦黏性系數來實現壁面函數,根據摩擦速度的定義,可以得到壁面剪切應力為

(5)

根據壁面剪切應力的定義：

(6)

式中:為壁面湍流黏性系數；和為鄰近壁面計算單元的切向速度和壁面距離,如圖2所示。根據式(5)、式(6)可以得到壁面湍流黏性系數為

(7)

基于國家數值風洞(NNW)軟件平臺,將SST-兩方程湍流模型和壁面函數模塊進行了耦合。NNW軟件平臺的基礎算法為有限體積方法,其中黏性項采用二階精度中心型格式離散,對流項采用二階TVD(Total Variation Diminishing)型格式,時間項采用LU-SGS隱式迭代算法。為了適應復雜外形,采用了分區網格算法。壁面邊界采用虛擬點方法實現,因此,通過修正虛擬點(如圖2中點)的值,來達到修正壁面剪切應力的目的。虛擬點湍流黏性系數為

(8)

以上標準壁面函數確定的虛擬點湍流黏性系數沒有涉及到具體的湍流模型,實際上,對邊界上湍流變量也需要進行適當的處理。壁面湍動能耗散率在黏性底層對數率層分別為

(9)

式中:渦黏性常數=0.09。

對于SST-兩方程湍流模型,壁面網格點上的湍流輸運變量可相應求出,其中

(10)

壁面虛擬點的湍流輸運變量為

(11)

求解SST-兩方程湍流模型的RANS方程時,結合以上描述的標準壁面函數,規定當>11.06時采用壁面函數邊界條件,當≤11.06時,SST-兩方程湍流模型直接作用到壁面。由于摩擦速度與壁面剪切應力相互耦合,利用壁面函數不能顯式求得計算等量,需要采用迭代方式,具體過程如圖3所示。

圖3 壁面函數程序實現流程圖Fig.3 Flow chart of wall function method in code

3 算例與分析

3.1 壓縮拐角

本算例用于考察壁面函數對湍流邊界層和分離流動的預測結果精度的影響規律。壓縮拐角來流馬赫數為2.85,來流壓力為23 600 N/m,來流溫度為100.3 K,等溫壁為276 K,湍流邊界層厚度=2.3 cm,壓縮拐角=24°。

圖4為計算采用的網格,為了能夠精確捕捉拐角處的激波邊界層干擾,對網格拐角處進行了加密,在拐角處網格流向距離為0.1 mm,整體規模為201×51(流向×法向),壁面距離分別取0.001 mm、0.01 mm、0.1 mm和1 mm,其中密網格=0.001 mm,分離區之前≈0.1,粗網格=1 mm,分離區之前≈200。

圖4 壓縮拐角網格示意圖Fig.4 Compression corner mesh

圖5給出了壁面距離=0.001 mm和1 mm網格計算消耗的計算時間和平均殘差RES之間的對比,從圖中可知粗網格=1 mm在下降到殘差10時,消耗的計算時間僅為密網格=0.001 mm 計算時間的1/4。

圖5 壓縮拐角24°收斂曲線Fig.5 Residual curves of 24° compression corner

圖6給出了不同網格間距湍流邊界層內速度分布,當不使用壁面函數時(如圖6(a)所示),隨著第1層壁面網格距離的增大,邊界層內速度分布和實驗值偏差逐步擴大,對于壁面網格間距=1 mm網格,和實驗值差異最大可以達到30%。圖6(b)給出了使用壁面函數后的速度分布,從圖中可知網格從最密的=0.001 mm逐漸變稀到=1 mm邊界層內速度分布依然和實驗值高度吻合,表明壁面函數的應用明顯提高了計算結果精度。

圖6 湍流邊界層速度分布Fig.6 Velocity distribution in turbulence boundary layer

圖7給出了不同網格間距時壓縮拐角壁面摩擦系數分布,從圖中可知=0.001 mm密網格分離前湍流邊界層摩擦系數同實驗值吻合較好,且可以準確模擬分離區起始位置,分離區再附位置稍落后于實驗值,再附區之后摩擦系數較實驗值有明顯偏差。文獻[27]使用雷諾應力模型(Reynolds Stress Model, RSM)和Launder-Sharma-湍流模型對該算例進行了數值模擬,再附區之后均出現了高估摩擦系數的現象,其中RSM模型預測值為實驗值的3倍左右。隨著網格間距增大,壁面摩擦系數同實驗值差異越大,其中粗網格=1 mm預測的摩擦系數在分離區前比實驗值低80%,嚴重低估摩擦系數。圖7(b)給出了使用壁面函數后壓縮拐角壁面摩擦系數分布,其中粗網格的預測能力相較圖7(a)得到大幅提升,分離區前后粗網格摩擦系數增加了10倍之多。

圖7 壓縮拐角摩擦系數分布Fig.7 Skin friction distribution for compression corner

3.2 高速飛行器

本算例用于檢驗標準壁面函數對于復雜的工程外形的氣動力預測能力。飛行器外形類似于美國X-51吸氣式高速飛行器,包含V尾、進氣道、機翼、安定面等部件。飛行馬赫數為4.0,雷諾數來流壓力2 634 N/m,來流溫度為69 K,絕熱壁,攻角范圍-8°～8°。為了保證數值預測精度,對進氣道唇口附近網格進行了十分精細的網格設計,進氣道內部采用O型網格。在飛行器表面生成若干層(包含41個點)較密的網格,安定面、尾舵、內流道等都采用了O型網格,在此基礎上再向空間拓展,半場密網格約900萬,第1層網格距離為2×10mm。稀網格將壁面O型網格修改為13個點,其他方向網格保持不變,第1層網格距離為0.1 mm,半場網格560萬,進氣道內密網格和稀網格的對比見圖8,其中圖中左邊界為對稱面、上邊界、下邊界和右邊界為壁面條件。

圖8 進氣道內網格對比Fig.8 Comparison of inlet meshes

圖9給出了高速飛行器對稱面馬赫數云圖和壓力等值線,展示進氣道內復雜的波系結構,且隨著攻角的增加,進氣道前緣產生的入射激波增強,入射激波和前體湍流邊界層干擾增強,在壁面附近產生的分離渦逐漸增大。

圖9 高速飛行器不同攻角馬赫數云圖和壓力等值線Fig.9 Pressure iso-lines superimposed on Mach number contour for high-speed flight vehicle with different angles of attack

附面層網格變稀后,僅對局部流場產生影響,其中入口處分離區大小發生明顯變化,圖10選取攻角=8°給出了高速飛行器唇口附近壓力云圖和流線。圖中:WF0表示不使用壁面函數;WF1表示使用壁面函數。從圖中可知唇口產生的入射激波和上表面的邊界層發生激波邊界層干擾,在上表面形成一個小的分離區,其中圖10(a)中給出的密網格結果可觀察到明顯的分離區,而稀網格不使用壁面函數時(見圖10(b))模擬出的分離區小得多,圖10(c)中稀網格使用壁面函數后得到的分離區大小和密網格更為接近,使得整個進氣道入口處的波系結構也更靠近密網格結果。

圖10 高速飛行器唇口附近壓力云圖和流線Fig.10 Stream lines superimposed on pressure coefficient contour for high-speed flight vehicle near inlet

圖11給出了進氣道內上表面壁面壓力系數對比。由于進氣道內存在復雜的波系結構,在進氣道上表面壓力系數呈現出多個波峰和波谷,稀網格不使用壁面函數時和密網格之間差距較大,表明稀網格不能準確捕捉進氣道內激波位置,而稀網格使用壁面函數后整體預測精度提高,其和密網格之間的偏差減小。圖中在進氣道入口處對壓力系數進行了放大,即圖10中相應位置,其中由于稀網格不使用壁面函數時不能預測壁面附近激波邊界層干擾產生的分離區,而使用壁面函數后不僅精確預測分離區位置,且可得到和密網格相似的壓力系數曲線,高壓區預測偏差從15%減少到2%。

圖11 進氣道內壓力系數曲線對比Fig.11 Pressure coefficient comparison in inlet

圖12給出了相同計算設置的情況下(僅附面層網格和壁面函數使用情況不同),計算相同的步數消耗的計算機CPU時間對比,從圖中可知減少附面層網格后計算時間減少了40%左右,而開啟壁面函數后計算量相對僅增加了1%。

圖12 相同計算設置消耗CPU時間Fig.12 CPU time consumption with the same numerical setup

圖13給出了軸向力隨CPU時間的收斂曲線,從圖中可知密網格整體軸向力的收斂較慢,軸向力系數處于緩慢下降狀態,稀網格整體收斂速度較快,整體看可大致節約60%的計算時間。

圖13 軸向力收斂曲線Fig.13 Convergence curves of axial force

圖14給出了軸向力系數隨攻角變化,隨著攻角增加軸向力系數增加,其中圖14(a)給出了軸向力系數摩擦力分量A-,摩擦力占整體比例為37%～50%,當粗網格不采用壁面函數時,嚴重低估摩擦力,和密網格結果相差40%左右,而使用壁面函數后,和密網格差異在4%之內。圖14(b)給出了軸向力系數壓力分量A-,從圖中可知不使用壁面函數時和密網格相差10%左右,使用壁面函數后相差3%,計算精度的提高主要得益于使用壁面函數后對分離區的模擬精度提高,如圖10 所示。整體軸向力如圖14(c)所示,整體軸向力粗網格不使用壁面函數時和密網格差異在11%左右,而使用壁面函數后差異縮小到2%以內,有效提高了軸向力預測精度。

圖14 高速飛行器軸向力隨攻角變化曲線Fig.14 Axial force variation with different angles of attack for high-speed flight vehicle

4 結 論

基于國家數值風洞可壓縮流動氣動力數值模擬平臺,開展適用的壁面函數方法研究,結合工程常用的SST-兩方程模型,通過子迭代的形式計算摩擦速度,在達到預定收斂準則后通過更新虛擬點湍流黏性系數來實現湍流模型與壁面函數的耦合,并通過典型二維算例進行驗證和工程復雜外形進行考核驗證,通過比較不同網格間距時壁面函數開關對數值結果的影響,得到結論如下：

1) 壁面函數可有效提高壁面速度型預測精度,從而提高壁面摩擦系數的預測精度,對于簡單外形摩擦系數大約可提高80%,對于工程復雜外形摩阻預測精度可提高36%。

2) 壁面函數可顯著放寬第1層壁面距離,減少網格總量,節約計算時間。對于簡單外形無量綱壁面距離≤200范圍內,均可準確模擬湍流邊界層速度分布。

3) 對于工程外形內外流一體化數值模擬,壁面函數在減少壁面網格距離和數量的同事,可準確模擬進氣道內復雜的激波邊界層干擾,激波反射等現象。

總的來說,通過壁面函數可獲得和密網格相近的預測結果同時,有效降低湍流模擬對第1層網格高度的限制,大幅節約計算時間和計算內存,但壁面函數在工程外形上的應用準則,還需進一步研究。由于隨著馬赫數的增加,Navier-Stokes方程中的動量方程和能量方程耦合作用加強,而本文的壁面函數可視為動量方程在壁面附近的近似解,未考慮能量方程的影響,因此在高馬赫數湍流模擬中,其預測精度很難得到保證,建議本文設計的壁面函數在來流馬赫數≤5.0范圍內的湍流模擬中使用。下一步將在國家數值風洞氣動力計算平臺中添加溫度壁面函數,提高粗網格氣動熱參數的預測精度,進一步擴展壁面函數的適用范圍。