壁面函數(shù)在超聲速湍流模擬中的應(yīng)用

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

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

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

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

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

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

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

本文基于國家數(shù)值風(fēng)洞氣動力計算平臺, 針對工程湍流模型SST-,開展適用的壁面函數(shù)研究,給出了具體的程序?qū)崿F(xiàn)過程和邊界處理形式,通過典型二維算例完成了壁面函數(shù)模塊考核,最后基于內(nèi)外流一體化的超聲速飛行器湍流流動的模擬,完成了模塊對復(fù)雜外流的工程應(yīng)用考核。

1 壁面函數(shù)

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

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

(1)

(2)

式中:

=()12=

(3)

=()12=

(4)

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

2 壁面函數(shù)在三維邊界層中的推廣應(yīng)用

從平面推廣到三維曲面,就是把邊界層方程及其相關(guān)理論表達式轉(zhuǎn)換到貼體坐標(biāo)系中,如圖2 所示。以當(dāng)?shù)仉x開壁面第1個計算單元處的切向速度作為流動方向,將三維流動進行局部二維化處理,在此基礎(chǔ)上采用湍流邊界層的壁面函數(shù),得到壁面剪切應(yīng)力。

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

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

(5)

根據(jù)壁面剪切應(yīng)力的定義：

(6)

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

(7)

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

(8)

以上標(biāo)準(zhǔn)壁面函數(shù)確定的虛擬點湍流黏性系數(shù)沒有涉及到具體的湍流模型,實際上,對邊界上湍流變量也需要進行適當(dāng)?shù)奶幚?。壁面湍動能耗散率在黏性底層對?shù)率層分別為

(9)

式中:渦黏性常數(shù)=0.09。

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

(10)

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

(11)

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

圖3 壁面函數(shù)程序?qū)崿F(xiàn)流程圖Fig.3 Flow chart of wall function method in code

3 算例與分析

3.1 壓縮拐角

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

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

圖4 壓縮拐角網(wǎng)格示意圖Fig.4 Compression corner mesh

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

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

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

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

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

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

3.2 高速飛行器

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

圖8 進氣道內(nèi)網(wǎng)格對比Fig.8 Comparison of inlet meshes

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

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

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

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

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

圖11 進氣道內(nèi)壓力系數(shù)曲線對比Fig.11 Pressure coefficient comparison in inlet

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

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

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

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

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

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

4 結(jié) 論

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

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

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

3) 對于工程外形內(nèi)外流一體化數(shù)值模擬,壁面函數(shù)在減少壁面網(wǎng)格距離和數(shù)量的同事,可準(zhǔn)確模擬進氣道內(nèi)復(fù)雜的激波邊界層干擾,激波反射等現(xiàn)象。

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