王 沛, 李 錦,2, 江定武,2, 毛枚良,2

(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) 621000)

引 言

對(duì)高超聲速飛行器而言, 隨著飛行高度增加, Mach數(shù)增加, Reynolds數(shù)減少, 表面邊界層厚度增加, 形成位移效應(yīng), 導(dǎo)致飛行器有效幾何外形發(fā)生變化, 影響邊界層外部的無(wú)黏流場(chǎng), 而外部流場(chǎng)的變化也會(huì)影響壁面附近的流場(chǎng)。邊界層與外部壓縮無(wú)黏流場(chǎng)之間這種強(qiáng)耦合作用形成黏性干擾效應(yīng)。這種效應(yīng)會(huì)引起當(dāng)?shù)貕毫Α?摩擦阻力以及熱傳導(dǎo)等參量的劇烈變化, 進(jìn)而影響飛行器的整體氣動(dòng)性能, 是高超聲速飛行時(shí)不可忽視的復(fù)雜物理效應(yīng)。

早期的黏性干擾研究主要采用理論分析的方法, 基于高超聲速邊界層相似律和高超聲速無(wú)黏流的壓力關(guān)聯(lián)公式[1]: Bertram[2]結(jié)合摩阻與壓力梯度的近似關(guān)系式, 研究了零厚度帶攻角平板和三角翼的黏性干擾效應(yīng)對(duì)壓力、 摩阻和氣動(dòng)力的影響; Whitfield等[3]完成了冷壁細(xì)長(zhǎng)鈍錐的黏性干擾效應(yīng)理論及與實(shí)驗(yàn)的對(duì)比研究; Stollery[4]針對(duì)指數(shù)律外形給出了凹陷體、 壓縮拐角和膨脹拐角等若干實(shí)例的黏性干擾流場(chǎng)結(jié)果; Probstein等[5]評(píng)估了鈍錐外形上橫向曲率在弱黏性干擾區(qū)的影響; Yasuhara[6]研究了3/4冪次律軸對(duì)稱外形上較弱至中等強(qiáng)度黏性干擾區(qū)中的相似流動(dòng)。但是早期這些研究工作基本上只針對(duì)簡(jiǎn)單外形, 復(fù)雜外形面臨更加復(fù)雜的流動(dòng)機(jī)理, 這些黏性理論可能失效, 難以得到正確的結(jié)果。

20世紀(jì)70年代, 美國(guó)對(duì)航天飛機(jī)OV102進(jìn)行了一系列飛行實(shí)驗(yàn), 獲得了大量的氣動(dòng)力飛行數(shù)據(jù)[7], 給出了與氣動(dòng)特性參量相關(guān)的第1、 第2、 第3黏性干擾參數(shù), 并得到第3黏性干擾參數(shù)數(shù)據(jù)關(guān)聯(lián)結(jié)果更為準(zhǔn)確的結(jié)論[8]?；诖? 國(guó)內(nèi)學(xué)者在黏性干擾效應(yīng)方面取得了一些研究成果。莊逢甘等[9]針對(duì)OV102外形進(jìn)行了初步研究, 提出了較為成熟的超聲速飛行器黏性干擾效應(yīng)的工程計(jì)算方法; 龔安龍等[10]對(duì)軌道器再入階段的黏性干擾效應(yīng)進(jìn)行了研究, 考察了黏性干擾相關(guān)參數(shù)的關(guān)聯(lián)特性; 毛枚良等[11]對(duì)升力體簡(jiǎn)化外形的高超聲速繞流流動(dòng)進(jìn)行了模擬, 認(rèn)為黏性干擾導(dǎo)致的氣動(dòng)力系數(shù)的改變量與黏性干擾參數(shù)在特定范圍內(nèi)具有良好的線性關(guān)系; 李維東等[12]結(jié)合參考溫度方法提出了一種能夠考慮黏性干擾效應(yīng)的高超聲速乘波體氣動(dòng)性能的工程預(yù)測(cè)方法; 陳堅(jiān)強(qiáng)等[13]針對(duì)OV102外形以及高升阻比復(fù)雜外形黏性干擾效應(yīng), 提出基于Mach數(shù)和第3黏性干擾參數(shù)的關(guān)聯(lián)參數(shù)和公式, 具有實(shí)現(xiàn)高超聲速氣動(dòng)力數(shù)據(jù)天地相關(guān)的能力; 張益榮等[14]針對(duì)典型高超聲速翼身組合體外形, 建立了完全氣體條件下縱向氣動(dòng)力系數(shù)的黏性干擾模型, 并進(jìn)行了不確定度量化分析; 龔安龍等[15]針對(duì)高升阻比面對(duì)稱飛行器, 建立了黏性干擾增量與黏性干擾參數(shù)之間的線性關(guān)系式, 并采用Mach數(shù)效應(yīng)、 真實(shí)氣體效應(yīng)等進(jìn)行修正建立了天地?fù)Q算公式; 羅長(zhǎng)童等[16]利用專業(yè)化機(jī)器學(xué)習(xí)算法, 得到了不同風(fēng)洞試驗(yàn)結(jié)果共同遵守的不變規(guī)律, 為建立黏性干擾等物理效應(yīng)天地?fù)Q算方法提供了新思路。

國(guó)內(nèi)外的研究現(xiàn)狀表明, 傳統(tǒng)黏性干擾的研究主要基于連續(xù)性假設(shè)成立的N-S方程, 但是真實(shí)飛行條件下隨著速度和高度進(jìn)一步變大, 氣體分子遠(yuǎn)離平衡態(tài), 稀薄效應(yīng)變得不可忽略, 導(dǎo)致飛行器的氣動(dòng)特性相比于連續(xù)流有了較大差異, 而這一區(qū)域內(nèi)目前仍然缺少合適的黏性干擾效應(yīng)預(yù)測(cè)模型, 這種模型的建立對(duì)求解器提出了新要求。統(tǒng)一氣體動(dòng)理論(unified gas kinetic scheme, UGKS)是近些年來(lái)比較熱門(mén)的一種跨流域多尺度求解方法, 在連續(xù)流區(qū)域可以得到精確的N-S解, 在稀薄區(qū)域也可以得到與DSMC一致的結(jié)果, 已經(jīng)廣泛應(yīng)用于從低速到高速、 從連續(xù)流到稀薄流的全速域全流域的流場(chǎng)模擬中[16]。

基于此, 本文采用UGKS求解器, 以典型外形鈍錐為研究對(duì)象。預(yù)測(cè)其在高空(70～110 km)、 高M(jìn)ach數(shù)(≥10)條件下的氣動(dòng)特性, 并借助黏性干擾效應(yīng)的理論成果, 將Mach數(shù)、 攻角、 黏性系數(shù)等參數(shù)與UGKS有黏解和Euler無(wú)黏解的差量進(jìn)行關(guān)聯(lián), 建立了氣動(dòng)特性增量預(yù)測(cè)模型。這樣后續(xù)進(jìn)行流場(chǎng)預(yù)測(cè)時(shí), 僅需要計(jì)算Euler無(wú)黏解, 就可通過(guò)該增量模型快速高效得到相應(yīng)黏性條件下的氣動(dòng)力參數(shù), 這極大簡(jiǎn)化了氣動(dòng)特性的求解過(guò)程。最后選取相關(guān)性曲線與Pearson乘積矩相關(guān)系數(shù)進(jìn)行了數(shù)據(jù)相關(guān)性分析, 采用部分新?tīng)顟B(tài)驗(yàn)證了模型預(yù)測(cè)結(jié)果能夠滿足工程應(yīng)用的精度要求。

1 數(shù)值方法

1.1 無(wú)黏解算器

控制方程選用一般曲線坐標(biāo)系下的三維可壓縮Euler方程, 基于有限體積法進(jìn)行離散, 采用隱式LUSGS方法進(jìn)行迭代求解, 具體計(jì)算過(guò)程參見(jiàn)文獻(xiàn)[17]。

1.2 黏性解算器

控制方程采用Shakhov模型方程, 其具體形式為

其中,f是氣體在空間x, 速度u, 以及時(shí)間t的分布函數(shù)。τ為碰撞時(shí)間, 與黏性和壓力有關(guān),f+可以表示為

其中, 第1項(xiàng)g是Maxwell分布函數(shù), 第2項(xiàng)是基于原始BGK模型方程的修正項(xiàng), 目的是為了得到合理的Prandtl數(shù),c和q分別為分子熱運(yùn)動(dòng)速度矢量和熱流矢量。

方程有如下積分解

e-t/τf(x-ut,0,u)

該積分解由初始分布函數(shù)和平衡態(tài)分布函數(shù)構(gòu)成, 通過(guò)松弛時(shí)間τ將分子的碰撞和自由輸運(yùn)過(guò)程耦合起來(lái), 從而實(shí)現(xiàn)分布函數(shù)在稀薄流和連續(xù)流區(qū)域的自動(dòng)調(diào)節(jié), 其具體表達(dá)式見(jiàn)文獻(xiàn)[18]。

宏觀守恒量Q以及應(yīng)力P、 熱流q等變量均可以通過(guò)在速度空間上對(duì)f求矩得到

其中, dΞ=dudvdw為相空間的體積元。

基于上述模型構(gòu)造了一套隱式三維大規(guī)模并行軟件, 并進(jìn)行了多個(gè)算例的驗(yàn)證, 詳見(jiàn)文獻(xiàn)[19-20]。

2 結(jié)果及建模

相較于返回艙, 鈍錐適用于需要大橫向機(jī)動(dòng)能力、 良好的機(jī)動(dòng)性、 低著陸誤差和低減速負(fù)載的任務(wù)以及高進(jìn)入速度和稀薄大氣層的任務(wù)。本文采用鈍錐外形進(jìn)行高超聲速黏性干擾模型的研究。鈍錐球頭半徑取為600 mm, 半錐角10°, 全長(zhǎng)L=3 600 mm, 力系數(shù)參考面積為1.0 m2, 力矩參考點(diǎn)位于質(zhì)心, 坐標(biāo)為(2 418 mm, 0 mm, 0 mm)。圖1給出了鈍錐的計(jì)算網(wǎng)格示意圖, 表1給出了計(jì)算狀態(tài)。

表1 計(jì)算狀態(tài)Table 1 Computation states

2.1 流場(chǎng)及氣動(dòng)力特性分析

圖2給出了Mach數(shù)10、 攻角10°時(shí)有黏(100 km)以及無(wú)黏兩個(gè)狀態(tài)的流場(chǎng)壓力云圖和速度矢量圖, 圖中變量均為無(wú)量綱結(jié)果。100 km時(shí), 黏性邊界層很厚, 黏性干擾效應(yīng)明顯改變了物面壓力和空間流場(chǎng)結(jié)構(gòu)。

(a) Inviscid flow

(b) Viscid flow(100 km)圖2 速度矢量及壓力云圖Fig. 2 Velocity vector and pressure contour

圖3給出了兩個(gè)高度下對(duì)稱面及物面局部Knudsen數(shù)分布。圖中局部Knudsen數(shù)變化范圍為0.01～100, 而且隨著高度進(jìn)一步增加, 激波內(nèi)部、 壁面附近的局部Knudsen數(shù)進(jìn)一步增加, 稀薄效應(yīng)更加明顯。這些局部區(qū)域內(nèi)傳統(tǒng)基于連續(xù)性假設(shè)的N-S方程已不再適用, 模擬的流場(chǎng)已不再準(zhǔn)確。為了更好地捕捉流場(chǎng)細(xì)節(jié), 預(yù)測(cè)真實(shí)流場(chǎng), 需要采用跨流域求解器。

(a) H=70 km

(b) H=80 km圖3 對(duì)稱面及物面局部Knudsen數(shù)分布Fig. 3 Local Kn number on symmetry plane and wall

圖4給出了Mach數(shù)為10、 攻角10°時(shí)不同高度對(duì)稱面固壁壓力分布比較, 在球頭與錐的交接位置(x/L=0.167)存在明顯的壓力分布拐點(diǎn)。迎風(fēng)面大面積區(qū)域(錐身)壓力分布隨著高度的增加呈增加趨勢(shì)。背風(fēng)面大面積區(qū)域壓力分布隨高度增加呈現(xiàn)先增后減的趨勢(shì), 但是量值較迎風(fēng)面大約要小1個(gè)量級(jí), 這與物理實(shí)際基本一致。

(a) Windward

(b) Leeawrd圖4 對(duì)稱面壓力比較Fig. 4 Comparison of pressure on symmetry plane

圖5, 6分別給出了Mach數(shù)為10、 攻角為10°和30°時(shí)壁面壓力和摩擦阻力對(duì)軸向力系數(shù)的貢獻(xiàn)隨黏性干擾參數(shù)的變化。隨著黏性干擾參數(shù)的增加, 軸向力系數(shù)增加, 壓力、 黏性的貢獻(xiàn)也在增加, 而且相比壓力項(xiàng), 黏性項(xiàng)的作用占主導(dǎo), 且對(duì)比發(fā)現(xiàn)攻角變化對(duì)軸向力的增量影響較小。

圖5 攻角10°時(shí)軸向力系數(shù)特性Fig. 5 Axial force coefficient characteristics at α=10°

圖6 攻角30°時(shí)軸向力系數(shù)特性Fig. 6 Axial force coefficient characteristics at α=30°

圖7給出了部分狀態(tài)的俯仰力矩特性。高度100 km(V′∞～0.33)及以下, 俯仰力矩變化緩慢, 符號(hào)為正, 壓心在質(zhì)心之前, 飛行器處于靜不穩(wěn)定狀態(tài)。高度進(jìn)一步增加到110 km, 攻角為10°， 20°情況下, 俯仰力矩符號(hào)變?yōu)樨?fù)數(shù), 壓心后移到質(zhì)心之后, 飛行器處于靜穩(wěn)定狀態(tài)。且Mach數(shù)對(duì)俯仰力矩的影響較小。

圖7 俯仰力矩系數(shù)特性Fig. 7 Pitching moment coefficient characteristics

2.2 模型建立

以往的研究成果表明, 第3黏性參數(shù)能夠較好地將來(lái)流條件與氣動(dòng)特性參量關(guān)聯(lián)起來(lái)[13], 所以本文基于第3黏性干擾參數(shù)建立黏性干擾預(yù)測(cè)模型。該參數(shù)表達(dá)式為

其中,T′表示邊界層內(nèi)的參考溫度。從式中可以看出, Mach數(shù)一定的情況下, 高度增加, Reynolds數(shù)減小, 黏性干擾參數(shù)增加。為了更直觀地表征黏性干擾參數(shù)與氣動(dòng)參量之間的關(guān)系, 所有的對(duì)比圖均采用黏性干擾參數(shù)為橫坐標(biāo)。

圖8給出了10°攻角對(duì)稱面上幾個(gè)典型流向位置處的壓力改變量的變化情況。迎風(fēng)面, 高度低于110 km時(shí), 壓力改變量與黏性干擾參數(shù)線性關(guān)系比較好。背風(fēng)面, 基本看不到較好的線性關(guān)系。這說(shuō)明, 即便是鈍錐這種簡(jiǎn)單的外形, 壓力改變量與黏性干擾參數(shù)之間也不能呈現(xiàn)線性關(guān)系。

(a) Windward

(b) Leeward圖8 壓力增量隨黏性干擾參數(shù)的變化Fig. 8 Variation of pressure increment with the viscous interaction parameter

以軸向力和俯仰力矩系數(shù)為例, 給出黏性干擾模型的分析建立過(guò)程?？梢灶愃频氐玫椒ㄏ蛄ο禂?shù)預(yù)測(cè)模型。

圖9給出了軸向力和俯仰力矩系數(shù)變化量隨第3黏性干擾參數(shù)的變化, 這里的系數(shù)變化量表示UGKS求解器得到的有黏解與Euler方程求解得到的無(wú)黏解的差量。

(a) Windward

(b) Leeward圖9 力系數(shù)增量隨黏性干擾參數(shù)的變化Fig. 9 Variation of force coefficient increment with the viscous interaction parameter

從圖中可以看到, 高度小于100 km時(shí), 軸向力系數(shù)變化量與第3黏性干擾參數(shù)呈較好的線性關(guān)系, 但是隨著高度的增加, 則逐漸偏離線性關(guān)系。因而, 本文采用二次多項(xiàng)式進(jìn)行擬合, 并參照黏性干擾模型的推導(dǎo)過(guò)程[13], 初步給出了力系數(shù)的黏性干擾模型

并進(jìn)一步假設(shè)

2.3 模型評(píng)估

采用45°相關(guān)線對(duì)結(jié)果進(jìn)行初步分析。圖10給出了黏性干擾模型預(yù)測(cè)數(shù)據(jù)(以model表示)與數(shù)值模擬計(jì)算數(shù)據(jù)(以UGKS表示)的相關(guān)性曲線, 不同攻角和Mach數(shù)條件下數(shù)據(jù)基本分布在相關(guān)線附近, 可見(jiàn)數(shù)據(jù)間相關(guān)性較好。

(a) ΔCA

(b) ΔCm圖10 黏性干擾模型與數(shù)值模擬預(yù)測(cè)結(jié)果相關(guān)性Fig. 10 Correlation between viscous interaction model and numerical simulation results

根據(jù)上述相關(guān)性曲線, 采用相對(duì)正交距離考察黏性干擾模型數(shù)據(jù)擬合精準(zhǔn)度

其中,di為數(shù)據(jù)點(diǎn)到相關(guān)性曲線的正交距離,xr為數(shù)據(jù)點(diǎn)在相關(guān)性曲線上投影對(duì)應(yīng)的橫坐標(biāo), 示意圖見(jiàn)圖11。圖中橫縱坐標(biāo)分別為干擾模型與UGKS模擬得到的力/力矩系數(shù)增量。軸向力和俯仰力矩系數(shù)的相對(duì)正交距離計(jì)算結(jié)果如圖12所示。對(duì)于軸向力系數(shù), 高度大于70 km時(shí), 模型預(yù)測(cè)的力系數(shù)變化量的相對(duì)偏差都低于3%; 70 km時(shí), 相對(duì)偏差最大可達(dá)14%, 但是此時(shí)變化量的量值(約為0.236)僅為基準(zhǔn)量(Euler無(wú)黏解, 約為1.18)的1/5。基準(zhǔn)量疊加變化量得到的軸向力系數(shù)預(yù)測(cè)值與UGKS計(jì)算值偏差大約僅為3%。俯仰力矩系數(shù)變化量的相對(duì)正交距離具有類似的情況, 在變化量較小的時(shí)候, 相對(duì)正交距離甚至可達(dá)28%。但是此時(shí)的變化量量值比基準(zhǔn)量量值(Euler無(wú)黏解)低1～2個(gè)量級(jí)。

圖11 正交距離示意圖Fig. 11 Illustration of orthogonal distance

(a) ΔCA

(b) ΔCm圖12 力系數(shù)相對(duì)正交距離Fig. 12 Relative orthogonal distance of the force coefficient

采用統(tǒng)計(jì)學(xué)中Pearson乘積矩相關(guān)系數(shù)進(jìn)行進(jìn)一步的相關(guān)性分析

最后, 為了初步評(píng)估模型的準(zhǔn)確性, 采用UGKS求解器和黏性干擾模型預(yù)測(cè)了Mach數(shù)20、 攻角10°, 高度分別為100 km和110 km時(shí)的氣動(dòng)力增量, 具體結(jié)果及相對(duì)偏差見(jiàn)表2, 本文建立的模型預(yù)測(cè)的軸向力增量與UGKS求解器模擬結(jié)果偏差小于1%, 俯仰力矩增量的預(yù)測(cè)誤差在10%以內(nèi)。初步說(shuō)明模型預(yù)測(cè)的氣動(dòng)增量與UGKS黏性求解器的結(jié)果基本一致, 采用黏性干擾模型可以較好地實(shí)現(xiàn)不同Mach數(shù)下的數(shù)據(jù)相關(guān), 可以為工程實(shí)踐提供一定的依據(jù)。

表2 模型預(yù)測(cè)與UGKS模擬結(jié)果對(duì)比Table 2 Comparison between model prediction and UGKS simulation

3 結(jié)論

本文基于統(tǒng)一氣體動(dòng)理學(xué)算法UGKS, 針對(duì)高超聲速流場(chǎng)下的鈍錐外形, 采用UGKS求解器得到黏性解, 采用Euler求解器得到無(wú)黏解, 兩者相減得到氣動(dòng)特性增量, 然后通過(guò)第3黏性干擾參數(shù)將不同高度、 Mach數(shù)和攻角下的氣動(dòng)特性增量進(jìn)行關(guān)聯(lián), 建立了氣動(dòng)力系數(shù)的黏性干擾模型, 初步得到以下結(jié)論:

(1)在全流域范圍內(nèi), 氣動(dòng)力系數(shù)增量與第3黏性干擾參數(shù)之間不能以簡(jiǎn)單的線性關(guān)系進(jìn)行描述;

(2)基于UGKS求解器建立的黏性干擾模型, 在高度較高時(shí)考慮了稀薄氣體效應(yīng)的影響, 相比N-S求解器預(yù)測(cè)結(jié)果更加真實(shí)準(zhǔn)確;

(3)利用第3黏性干擾參數(shù)和來(lái)流參數(shù), 結(jié)合Euler無(wú)黏解以及全流域的黏性解, 可以建立具有良好精度的黏性干擾模型, 有助于高效獲得飛行器全流域氣動(dòng)力特性。

本文研究的鈍錐飛行器外形比較簡(jiǎn)單, 與真實(shí)飛行器有一定差別。把本文的建模方法擴(kuò)展至外形更加復(fù)雜的飛行器上是下一步努力的方向。

致謝審稿人提出的寶貴意見(jiàn)對(duì)提高文章質(zhì)量及完整性起到了很大的幫助和促進(jìn)作用, 在此表示感謝。