分子黏性對解析湍流壁面函數(shù)的影響

王新光, 陳琦, 萬釗, 高曉成, 燕振國

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

0 引言

高超聲速湍流邊界層模擬時,如需精確模擬壁面摩阻和熱流,要求壁面網格足夠密,通常無量綱壁面距離y+≈1,導致壁面附近網格驟增,使得計算收斂速度變慢,同時還會影響數(shù)值計算穩(wěn)定性。因此,主流商業(yè)軟件例如Fluent和CFD++,對工程外形進行湍流模擬時通常使用壁面函數(shù),大幅放寬壁面網格尺度,以提高計算效率。

壁面函數(shù)最初基于Prandtl混合長度理論和若干假設,得到不可壓縮湍流邊界層內的相似解[1],在湍流計算中通過引入壁面函數(shù)達到提高粗網格氣動力預測精度的目的[2]。在標準壁面函數(shù)的基礎上,文獻[3]通過對摩擦速度構造局部迭代實現(xiàn)壁面函數(shù)在超聲速流動中的應用,其中在復雜外形氣動力預測中可有效提高粗網格的預測精度。此外,通過添加壓力梯度[4]、旋轉修正[5]、可壓縮性和壁面?zhèn)鳠嵝猍5]等影響,進一步擴展了壁面函數(shù)在可壓縮湍流邊界層中的應用。文獻[6]通過使用Crocco-Busemann溫度方程,添加可壓縮性對壁面函數(shù)的影響,并將速度和溫度壁面函數(shù)耦合求解,發(fā)展了適用于可壓縮湍流邊界層的壁面函數(shù)。國內學者[7-8]將這種壁面函數(shù)耦合到了k-ω(k表示湍動能,ω為比耗散率)兩方程模型中,對于無量綱壁面距離y+<100范圍,取得滿意的結果。文獻[9]比較了兩方程k-ωSST、線性Launder-Sharmak-ε(ε表示湍動能耗散率)模型(LS)k-ε和非線性k-ε模型在超聲速和高超聲速流動中的應用,比較發(fā)現(xiàn)上述k-ω和k-ε模型的結果是類似的。文獻[10-11]基于文獻[6]中的壁面函數(shù),使用數(shù)值實驗和風洞實驗數(shù)據(jù)開展了修正研究,通過對速度、溫度的壁面函數(shù)修正來提高其預測精度。目前大多壁面函數(shù)研究的適用性僅針對湍流邊界層,對于存在分離和再附等逆壓梯度的復雜流動,適用性仍無法確定[12-14]。

近年來發(fā)展的解析壁面函數(shù)[15]對于存在分離的流動表現(xiàn)良好[16-18]。文獻[19]考慮壁面網格內對流項變化和能量方程中黏性耗散項的影響,發(fā)展了適用于可壓縮流動解析壁面函數(shù)(MAWF),通過二維超聲速激波邊界層干擾算例進行驗證,數(shù)值結果表明發(fā)展的可壓縮修正解析壁面函數(shù)消除了原始解析壁面函數(shù)的非物理振蕩,且大幅提升了壁面函數(shù)壁面熱流的預測精度,接近密網格低雷諾數(shù)模型結果,且可節(jié)約大量計算時間,文中馬赫數(shù)為5的算例其計算時間僅為密網格的5%。

可壓縮湍流邊界層直接數(shù)值模擬(DNS)結果[20]和標準壁面函數(shù)相比,當來流馬赫數(shù)較小時標準壁面函數(shù)依然適用,與DNS結果吻合,但隨著馬赫數(shù)的逐漸增大,標準壁面函數(shù)和DNS之間的偏差逐漸增大,表明低雷諾數(shù)超聲速槽道湍流存在顯著的可壓縮性。對于Ma為2.48的槽道流動,在黏性底層邊緣附近y+≈10,相較于壁面參數(shù),密度減小50%,溫度增高60%。通常黏性底層內溫度的快速變化必然導致流體黏性系數(shù)的變化,這一點對于壁面函數(shù)尤為重要。

本文基于發(fā)展的可壓縮解析壁面函數(shù),通過構造不同的黏性系數(shù)方程,考慮密度、黏性系數(shù)等可壓縮流動參數(shù)對無量綱壁面距離的影響,構造了兩種不同黏性系數(shù)分布,發(fā)展了兩種適應于高超聲速流動的解析壁面函數(shù),以達到提高解析壁面函數(shù)粗網格壁面熱流預測精度的研究目的,并通過高超聲速激波邊界層干擾算例,比較不同黏性系數(shù)方程對數(shù)值模擬結果的影響,為解析壁面函數(shù)的工程化應用起到拋磚引玉的作用。

1 可壓縮流動無量綱壁面距離

解析壁面函數(shù)[15,19]理論中,無量綱壁面距離的定義對于解析壁面函數(shù)至關重要。原始的解析壁面函數(shù)中無量綱壁面距離的定義為

(1)

式中:ρw表示壁面密度;y為壁面距離;kP為湍動能,下標P表示壁面第1層網格點;μw表示壁面黏性系數(shù)。

隨著馬赫數(shù)的增加,黏性底層內溫度、密度等參數(shù)發(fā)生很大變化,文獻[20]中通過對可壓縮湍流邊界層的DNS結果分析,認為可壓縮流動中無量綱壁面距離采用當?shù)刂迪噍^于壁面值更合適。但是壁面函數(shù)使用粗網格,壁面第1層網格內通常包含黏性底層和全湍流層,因此一種方式是使用主網格計算得到的當?shù)豍點的值定義無量綱壁面距離:

(2)

式中:ρP表示壁面第1層網格點密度;μP表示壁面第1層網格點黏性系數(shù)。

數(shù)值實驗結果表明,這種定義方式在計算激波邊界層干擾算例時穩(wěn)定性較差,為避免類似的問題出現(xiàn),本文采用黏性底層邊緣處的值定義,即

(3)

式中:ρv表示黏性底層位置的邊緣密度;μv表示黏性底層位置的黏性系數(shù)。

黏性底層處的密度和黏性系數(shù)使用上一時間步的解析溫度Tv計算,即后文式(9)在黏性底層邊緣處的值計算。使用這種方式定義的無量綱壁面距離避免了數(shù)值不穩(wěn)定性,同時考慮了壁面網格內流動參數(shù)變化的影響。

2 分子黏性系數(shù)建模

隨著來流馬赫數(shù)的增加,可壓縮湍流邊界層內溫度梯度增加(見圖1),圖1中δ為邊界層厚度,T為溫度,Tw為壁面溫度,Ma為來流馬赫數(shù),T∞為來流溫度,湍流邊界層來流條件可參考文獻[21-23]。而發(fā)展的可壓縮壁面函數(shù)[19],僅考慮壁面網格內湍流黏性系數(shù)的變化,可能導致解析壁面函數(shù)對于高超聲速流動的預測精度降低。

圖1 可壓縮平板湍流邊界層溫度分布

對于完全氣體,黏性系數(shù)僅是溫度的函數(shù),但直接將溫度解析表達式和黏性系數(shù)通過Sutherland公式連接后,很難得到解析的速度和溫度表達式。本文通過解析壁面函數(shù)得到黏性底層的溫度Tv,通過Sutherland公式得到黏性底層位置的黏性系數(shù)μv,并與壁面黏性系數(shù)關聯(lián),構造黏性底層內的黏性系數(shù)表達式。其中一種簡單的思路是與壁面處μw線性連接,如圖2(a)所示。圖2中,N為壁面第2個點,P為壁面第1個點,S為壁面中心,μN為N點的黏性系數(shù),yv為黏性底層的距離。

圖2 分子黏性系數(shù)分布示意圖

則黏性底層內的黏性系數(shù)為

(4)

當使用式(4)代入解析壁面函數(shù)中時,會使得解析溫度的表達式更加復雜,并使得解析壁面函數(shù)的穩(wěn)定性變差[24]。考慮到上述問題,文獻[24]對不同類型的黏性系數(shù)表達式進行了探索研究,其中兩種方式的黏性系數(shù)表達式在數(shù)學上接近于線性表達式(見式(1)),但不存在線性表達式的穩(wěn)定性問題,如圖2(b)和圖2(c)所示。

圖2(b)表示雙曲型分布:

(5)

圖2(c)表示拋物型分布:

(6)

在全湍流區(qū)使用μv,如圖2所示。在實際程序中僅需保存上一時間步的解析溫度Tv。

圖3給出了用式(5)和式(6)構造的解析壁面函數(shù)模擬Ma=8.18的激波邊界層算例(來流條件見表1,其中β為激波產生器角度,θ0為動量邊界層厚度,p∞為來流壓力)。干擾區(qū)前后壁面網格內溫度和分子黏性系數(shù)分布,并與密網格LS結果進行了對比,其中粗網格第1層網格高度為1.2 mm。由圖3可知,激波邊界層干擾區(qū)后黏性底層法向距離變小,且黏性底層內溫度變化梯度更大,而全湍流區(qū)的溫度變化較小,對應的分子黏性系數(shù),在黏性底層邊緣大約較壁面值增加了30%,其中采用式(5)的雙曲型分布時,分子黏性系數(shù)呈上凸形式,采用式(6)的拋物型分布時,呈下凹形式,且更接近于密網格低雷諾數(shù)模型結果,而壁面網格內兩種分布形式解析溫度分布差異很小,幾乎可以忽略,溫度分布結果接近于密網格結果,而文獻[19]中發(fā)展的MAWF其溫度分布和密網格有明顯差異。

表1 Ma=8.18時2D斜激波邊界層干擾來流條件

圖3 Ma=8.18算例粗網格解析溫度(左)和分子黏性系數(shù)(右)分布與密網格低雷諾數(shù)模型結果對比

3 雙曲型分布可壓縮解析壁面函數(shù)

解析壁面函數(shù)通過對邊界層內簡化的動量方程和能量方程進行積分,分別得到邊界層內速度和溫度表達式[17-19],當黏性系數(shù)采用雙曲型分布(見式(5))時,通過對簡化的動量方程進行兩次積分,可得到速度的解析表達式如下:

(7)

(8)

式中:常數(shù)α=0.229 5;下標1表示黏性底層(即y

vp為第1層網格點y方向的速度,p表示壓力;系數(shù)A1、A2、B2分別為

注意到本文發(fā)展的解析壁面函數(shù)公式中bμ=0時,將不考慮黏性系數(shù)的分布,使用黏性底層的μv值來代替壁面值μw,其公式與可壓縮解析壁面函數(shù)[18]相同。

對簡化的能量方程進行類似的積分,可得到解析溫度表達式為

(9)

(10)

式中:Pr為普朗特數(shù);cp為定壓比熱,cp=1 004.06 J/(kg·K);αt為常數(shù)0.193 6;Dth1、Dth2為能量方程的對流項系數(shù),

注意到解析溫度表達式(式(9)和式(10))和解析速度表達式(式(7)和式(8))耦合在一起,因此在程序中使用數(shù)值積分計算。這里僅給出雙曲型分布的最終表達式,拋物型的最終解析表達式可參考文獻[25]。

4 數(shù)值結果與分析

本文使用不同來流馬赫數(shù)的激波/湍流邊界層干擾算例來驗證考慮可壓縮流體邊界層參數(shù)變化的解析壁面函數(shù)、MAWF,與本文發(fā)展的雙曲型(hyper-MAWF)或拋物型(para-MAWF)分布解析壁面函數(shù)進行對比。

數(shù)值計算中無黏通量使用Roe-Pike格式,黏性通量采用2階中心差分格式。粗網格第1層網格y+大約為30,湍流模型使用標準k-ε模型,壁面模型采用上文描述的3種壁面函數(shù)。密網格第1層網格的y+≈1,湍流模型采用添加Yap修正[26]的LSk-ε模型。來流入口給定充分發(fā)展的湍流邊界層,具有1.5%的湍流度和黏性系數(shù)比率μT/μ=10。本文算例網格無關性研究表明計算結果與網格無關,詳情可參考文獻[25],其中Ma=8.18算例密網格和粗網格網格無關性結果如圖4所示,其中Q表示熱流,Qw表示平板壁面熱流。

圖4 Ma=8.18算例網格無關性研究

4.1 Ma=8.18斜激波邊界層干擾

斜激波邊界層干擾作為超聲速飛行器唇口處的簡化模型,斜激波以不同角度入射平板,會產生激波反射、分離和再附等復雜流動現(xiàn)象。Ma=8.18斜激波邊界層干擾[23],在激波邊界層干擾區(qū)域前局部來流條件如表1所示。由于激波產生器產生的斜激波和拐角的膨脹波對壁面均有干擾,在計算中上邊界給定激波產生器的壁面邊界模擬入射激波和激波產生器拐角的膨脹波,計算網格如圖5所示。

圖5 Ma=8.18算例網格示意圖

圖6給出了密網格480×210使用LS模型與粗網格160×60壁面函數(shù)方法湍流邊界層內的速度和溫度分布,可見數(shù)值解和實驗值基本吻合,其中密網格結果和實驗值吻合較好,在湍流邊界層內無量綱溫度呈現(xiàn)先減小后增大的趨勢,粗網格在邊界層內網格點較少,僅有10個點,大致可以刻畫出邊界層內的溫度變化。

圖6 Ma=8.18算例邊界層厚度0.94 mm位置邊界層內速度和溫度分布

圖7給出了β=10°時算例不同位置近壁面解析溫度與密網格LS模型的對比,其中x分別取值25 cm和55 cm時,截面位置位于激波邊界層干擾區(qū)前后,為全湍流邊界層。從圖7中可知,MAWF溫度梯度明顯較大,表明MAWF計算得到的壁面熱流值相對較大,在圖9壁面熱流分布中可對照觀察到。當x分別取值35 cm和45 cm時,即激波邊界層干擾區(qū)內,此時近壁面溫度梯度較大,本文設計的兩種壁面函數(shù)即hyper-MAWF和para-MAWF,溫度梯度相對其他兩種方法較小,對應圖9壁面熱流在分離區(qū)內較小。

圖7 Ma=8.18,β=10°算例解析溫度分布

圖8給出了β=10°算例不同位置近壁面解析速度分布,其中在激波干擾前后位置x=25 cm和x=55 cm處,完全湍流區(qū)處速度分布有較明顯的差異,其中MAWF解析速度梯度明顯較其他數(shù)值結果大,雖然圖6中邊界層4種數(shù)值方法計算得到的速度分布大致與實驗相同,但由于解析速度在壁面處的梯度決定了摩阻的數(shù)值,即MAWF計算得到更大的摩阻,而本文設計的兩種壁面函數(shù)其解析速度分布幾乎重合在一起,整體速度梯度較LS和MAWF較小。

圖8 Ma=8.18,β=10°算例解析速度分布

圖9比較了Ma=8.18算例不同激波產生器壁面壓力、摩擦阻力Cf和熱流結果(其中Q∞表示來流熱流),相較于MAWF,考慮流體性質的解析壁面函數(shù)得到的壁面壓力在分離區(qū)略有下降,壁面熱流在分離區(qū)顯著降低,更接近于實驗值,兩種不同分布形式黏性系數(shù)對應的壁面函數(shù)也出現(xiàn)差異,其中使用雙曲型分布得到的結果更低一些,相較于拋物型分布低約4%,而相較于MAWF降低了約35%。密網格LS結果出現(xiàn)高估壁面熱流現(xiàn)象,與文獻[27]中密網格低雷諾數(shù)模型結論一致。壁面摩阻和壁面熱流有相似的分布,除MAWF外,其他數(shù)值計算出現(xiàn)一個較小的分離區(qū)(摩阻為負),且在整個干擾區(qū)內本文設計的兩種壁面函數(shù)計算得到的最大摩阻比LS結果小約30%,有較大差異,但本文實驗沒有對應的摩阻試驗結果,暫時還無法確定壁面摩阻數(shù)值計算的精度。

圖9 Ma=8.18算例激波產生器5°壁面壓力(上)、摩阻(中)和熱流(下)分布

4.2 Ma=9.22壓縮拐角

壓縮拐角流動常見于高超聲速飛行器控制舵面、翼面和飛行器表面之間,本文選取Ma=9.22壓縮拐角[28],計算網格較為簡單,入口給定充分發(fā)展的湍流邊界層,來流參數(shù)如表2所示,其中Re表示雷諾數(shù)。其中低雷諾數(shù)模型使用的密網格為240×150,壁面函數(shù)使用的粗網格為180×80。其中密網格壓縮拐角角度β分別取值15°和32°時的馬赫數(shù)云圖如圖10所示,從中可知隨著馬赫數(shù)的增大,壓縮拐角產生的斜激波逐漸增強,與來流湍流邊界層之間的干擾加劇。

表2 Ma=9.22壓縮拐角來流條件

圖10 Ma=9.22算例馬赫數(shù)云圖

圖11比較了壓縮拐角分別為15°、32°、34°和38°時不同模型的壁面熱流分布,其中橫坐標S表示距離拐角處的位置(cm),低雷諾數(shù)LSk-ε模型高估了壁面熱流,當壓縮拐角38°時最大熱流幾乎是實驗值的3倍。文獻[29]使用LSk-ε模型和Wilcoxk-ω模型,文獻[30]使用Rodik-ε模型模擬了34°壓縮拐角。這3種低雷諾數(shù)模型均嚴重高估了壁面熱流值。MAWF也返回較高熱流值,在壓縮拐角上接近于LS模型結果,均高于實驗值,而本文發(fā)展的兩種壁面函數(shù)均給出合理的壁面熱流值,更接近于試驗結果,兩種黏性系數(shù)數(shù)值結果差異較小,在5%之內,且雙曲型分布更接近于實驗值。

圖11 Ma=9.22壓縮拐角壁面熱流

5 結論

本文基于高超聲速湍流邊界層特征,使用解析溫度計算黏性底層邊緣密度和黏性系數(shù),重新定義了無量綱壁面距離,并構造了雙曲型和拋物型兩種黏性系數(shù)方程,將黏性底層內快速的速度梯度變化考慮到解析壁面函數(shù)中,構造了拋物型解析壁面函數(shù)para-MAWF和雙曲型解析壁面函數(shù)hyper-MAWF,通過高超聲速邊界層干擾算例對構造的壁面函數(shù)進行了測試,并與低雷諾數(shù)LS和MAWF數(shù)值結果進行了對比。得出主要結論如下:

1)對于高超聲速算例,本文構造的解析壁面函數(shù)預測的壁面熱流結果更接近試驗結果,密網格LS模型和MAWF數(shù)值結果在激波邊界層干擾區(qū)均出現(xiàn)高估熱流現(xiàn)象。

2)本文構造的兩種不同黏性系數(shù)表達式,數(shù)值模擬結果差異較小,在5%之內,其中雙曲型分布預測的壁面熱流更接近于實驗值。

3)本文構造的解析壁面函數(shù)預測的速度、溫度和黏性系數(shù)分布更接近于密網格結果,而MAWF未考慮黏性系數(shù)分布,因此黏性系數(shù)和密網格之間差異較大。

整體來看,本文構造的解析壁面函數(shù)相較于MAWF壁面熱流的預測結果更加精確,在干擾區(qū)甚至顯示出優(yōu)于密網格LS模型的預測結果,且壁面函數(shù)采用粗網格進行數(shù)值模擬,可顯著節(jié)約計算時間[19]。由于兩種壁面函數(shù)數(shù)值計算結果差異較小,且形式上雙曲型更為簡單,對于高超聲速算例更推薦使用雙曲型分布的解析壁面函數(shù)。