王新光,毛枚良,何琨,*,陳琦,萬(wàn)釗

1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所,綿陽(yáng) 621000 2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,綿陽(yáng) 621010

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

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

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

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

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

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

1 壁面函數(shù)

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

黏性底層是一個(gè)緊貼固體壁面的極薄層,分子黏性切應(yīng)力是主要特征因素,湍流切應(yīng)力極小且可以忽略,所以流動(dòng)似乎處于層流狀態(tài),沿壁面法向呈線性分布,其區(qū)域?yàn)?。對(duì)數(shù)率層是一個(gè)距壁面較遠(yuǎn)的流體層,流動(dòng)處于完全湍流狀態(tài),湍流渦黏性遠(yuǎn)大于分子黏性,速度分布接近于對(duì)數(shù)率,位于從>到0.2倍邊界層厚度范圍內(nèi),在局部平衡的湍流邊界層內(nèi),可通過(guò)理論推導(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)理論表達(dá)式轉(zhuǎn)換到貼體坐標(biāo)系中,如圖2 所示。以當(dāng)?shù)仉x開(kāi)壁面第1個(gè)計(jì)算單元處的切向速度作為流動(dòng)方向,將三維流動(dòng)進(jìn)行局部二維化處理,在此基礎(chǔ)上采用湍流邊界層的壁面函數(shù),得到壁面剪切應(yīng)力。

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

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

(5)

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

(6)

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

(7)

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

(8)

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

(9)

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

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

(10)

壁面虛擬點(diǎn)的湍流輸運(yùn)變量為

(11)

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

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

3 算例與分析

3.1 壓縮拐角

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

圖4為計(jì)算采用的網(wǎng)格,為了能夠精確捕捉拐角處的激波邊界層干擾,對(duì)網(wǎng)格拐角處進(jìn)行了加密,在拐角處網(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)格計(jì)算消耗的計(jì)算時(shí)間和平均殘差RES之間的對(duì)比,從圖中可知粗網(wǎng)格=1 mm在下降到殘差10時(shí),消耗的計(jì)算時(shí)間僅為密網(wǎng)格=0.001 mm 計(jì)算時(shí)間的1/4。

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

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

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

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

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

3.2 高速飛行器

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

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

圖9給出了高速飛行器對(duì)稱(chēng)面馬赫數(shù)云圖和壓力等值線,展示進(jìn)氣道內(nèi)復(fù)雜的波系結(jié)構(gòu),且隨著攻角的增加,進(jìn)氣道前緣產(chǎn)生的入射激波增強(qiáng),入射激波和前體湍流邊界層干擾增強(qiáng),在壁面附近產(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)格變稀后,僅對(duì)局部流場(chǎng)產(chǎn)生影響,其中入口處分離區(qū)大小發(fā)生明顯變化,圖10選取攻角=8°給出了高速飛行器唇口附近壓力云圖和流線。圖中:WF0表示不使用壁面函數(shù);WF1表示使用壁面函數(shù)。從圖中可知唇口產(chǎn)生的入射激波和上表面的邊界層發(fā)生激波邊界層干擾,在上表面形成一個(gè)小的分離區(qū),其中圖10(a)中給出的密網(wǎng)格結(jié)果可觀察到明顯的分離區(qū),而稀網(wǎng)格不使用壁面函數(shù)時(shí)(見(jiàn)圖10(b))模擬出的分離區(qū)小得多,圖10(c)中稀網(wǎng)格使用壁面函數(shù)后得到的分離區(qū)大小和密網(wǎng)格更為接近,使得整個(gè)進(jìn)氣道入口處的波系結(jié)構(gòu)也更靠近密網(wǎng)格結(jié)果。

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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