彭傲平 李志輝 吳俊林 皮興才 蔣新宇

* (中國空氣動力研究與發(fā)展中心超高速空氣動力研究所,四川綿陽 621000)

? (北京航空航天大學(xué)國家計算流體力學(xué)實驗室,北京 100191)

引言

近年來,隨著對臨近空間戰(zhàn)略意義認識的加深,各航空航天大國發(fā)展了眾多臨近空間飛行器.為滿足長航時、大載荷、高超聲速巡航等需求,這類臨近空間飛行器通常整體尺寸較大,各部件間差異顯著,飛行過程會出現(xiàn)多物理多場耦合、多流區(qū)共存的多尺度復(fù)雜流動現(xiàn)象,傳統(tǒng)的數(shù)值模擬方法如基于宏觀連續(xù)介質(zhì)假設(shè)的納維?斯托克斯(N-S) 方程、模擬微觀分子運動與碰撞的直接模擬蒙特卡洛(DSMC)方法難以滿足多尺度一體化模擬需求,需要研究發(fā)展新的模擬手段[1-4].本文從介觀氣體分子速度分布函數(shù)滿足的方程即玻爾茲曼方程[5-7]出發(fā)開展多尺度流動一體化模擬研究.

玻爾茲曼方程作為氣體動理學(xué)理論的基本方程,能描述各流域氣體流動輸運現(xiàn)象,但它是一個高度非線性的積分?微分方程,其理論解和直接數(shù)值求解均難以實現(xiàn)[7].為利用玻爾茲曼方程能統(tǒng)一描述各流域氣體流動的特性,眾多學(xué)者提出了系列替代求解技術(shù),其中應(yīng)用最廣泛的是模型方程方法,即用一個相對簡單的表達式替換玻爾茲曼方程碰撞積分項[8],再對簡化后的偏微分方程進行求解.

模型方程作為玻爾茲曼方程的近似,必須滿足玻爾茲曼方程的基本屬性,如總和不變量性質(zhì)、H 定理、正確的輸運系數(shù)等.第一個也是最著名的模型方程是巴特納?格羅斯?克魯克模型[8](BGK模型),但通過對其進行查普曼?恩斯科(Chapman?Enskog)展開分析得到的普朗特數(shù)(Pr)為1,而單原子氣體的Pr約為2/3[9].為克服這一缺陷,國內(nèi)外學(xué)者提出了一系列模型方程,包括橢球統(tǒng)計?BGK (ESBGK)模型[10-11]、沙克霍夫(Shakhov)模型[12-13]、劉(Liu)模型[14]、貝利(Belyi)模型[15]等,這些模型方程均假定其碰撞頻率與分子速度無關(guān).此外,也有學(xué)者發(fā)展了考慮碰撞頻率與分子速度相關(guān)性的模型方程[16]以及適于多原子氣體、考慮熱力學(xué)非平衡的模型方程[17].

大量數(shù)學(xué)分析表明[16,18-20],上述模型方程均能滿足玻爾茲曼方程的基本屬性,且與原始玻爾茲曼方程的碰撞項具有類似的松弛特性,能保持與納維?斯托克斯方程、直接模擬蒙特卡洛方法的收斂一致,因結(jié)構(gòu)簡單、貼近跨流域流體物理變化特點、可計算技術(shù)針對性強而得到廣泛應(yīng)用[21],如泰塔雷夫(Titarev)建立了物理空間和速度空間均采用非結(jié)構(gòu)網(wǎng)格的內(nèi)斯韋塔依(Nesvetay)求解器[22-23],針對返回式空間飛行器等高馬赫數(shù)流動與直接模擬蒙特卡洛方法進行了系列對比分析.李志輝等[24-31]從研究描述各流域均適用的氣體分子速度分布函數(shù)方程出發(fā),吸收計算數(shù)學(xué)指數(shù)型積分求解原理,發(fā)展了氣體分子運動論離散速度坐標法,研制經(jīng)改進的高斯?埃爾米特積分方法等系列離散速度數(shù)值積分技術(shù),消除氣體分子速度分布函數(shù)對速度空間的連續(xù)依賴性,將玻爾茲曼模型方程化為在各個離散速度坐標點處具有非線性源項的雙曲型守恒方程,基于計算流體力學(xué)數(shù)值求解技術(shù),構(gòu)造直接求解速度分布函數(shù)演化更新的氣體動理論數(shù)值格式和大規(guī)模并行計算策略,建立了氣體動理論統(tǒng)一算法(GKUA),開展了自由分子流到連續(xù)流復(fù)雜高超聲速氣動力/熱問題大規(guī)模并行計算及其在航天工程中的應(yīng)用研究,同時對跨流域多體干擾流動問題、跨流域非定常流動問題進行了研究,獲得較好的成功.徐昆等[32-37]利用上述離散速度坐標法與BGK 類型格式相耦合,發(fā)展了統(tǒng)一氣體動理學(xué)格式(UGKS),使之適于模擬稀薄流到連續(xù)流氣體流動問題.

另一方面,從流動反映的物理現(xiàn)象來看,當流場中出現(xiàn)非平衡現(xiàn)象時,某點的速度分布函數(shù)已經(jīng)不再是簡單的麥克斯韋平衡態(tài)分布,而是呈現(xiàn)各種形態(tài),如雙峰分布、偏心分布等[38].這些分布與當?shù)氐臍怏w分子速度、流場宏觀速度、溫度以及熱流矢量、應(yīng)力張量等宏觀輸運性質(zhì)相關(guān),這從玻爾茲曼方程的查普曼?恩斯科分析解[5,39]可以看出.因此在對玻爾茲曼方程碰撞積分進行物理分析、構(gòu)造碰撞項的模型方程時應(yīng)同時考慮這些因素,而常用的玻爾茲曼碰撞積分模型如沙克霍夫模型沒有考慮應(yīng)力張量的影響,ES-BGK 模型屬于各向異性的類麥克斯韋分布,沒有考慮熱流矢量的影響.當模型方程不能全面反映這些因素時會導(dǎo)致模擬結(jié)果出現(xiàn)系統(tǒng)誤差,例如在稀薄過渡流區(qū)流動的克努森層內(nèi)黏性效應(yīng)占主導(dǎo),應(yīng)考慮應(yīng)力張量的影響,沙克霍夫模型與直接模擬蒙特卡洛方法模擬結(jié)果存在較大差別[40].

為此,本文將從玻爾茲曼方程出發(fā),基于查普曼?恩斯科漸近分析解[9]的數(shù)學(xué)推導(dǎo),構(gòu)造一種同時考慮熱流矢量和應(yīng)力張量影響、滿足玻爾茲曼方程高階碰撞矩的綜合型可計算碰撞模型,并在數(shù)學(xué)上分析其基本屬性,建立與玻爾茲曼方程的聯(lián)系,給出新模型與現(xiàn)有模型方程的關(guān)系,同時基于碰撞動力學(xué)確定該可計算模型碰撞松弛參數(shù)表達式.最后在氣體動理論統(tǒng)一算法框架下,使用該新建玻爾茲曼方程可計算模型,通過模擬一維激波結(jié)構(gòu)、二維近空間飛行環(huán)境平板和多體圓柱干擾流動來驗證新模型的有效性和可靠性.

1 模型方程的構(gòu)造

其中,i=x,y,z,Vx,Vy,Vz為氣體分子速度分量;x,y,z為空間位置坐標,t為時間;S(f),(?f/?t)col為碰撞積分項.

在確定速度分布函數(shù)f后,可以通過對其取矩得到宏觀流動參數(shù),如密度 ρ,宏觀速度分量Ux,Uy,Uz,溫度T,應(yīng)力張量分量Pij;熱流矢量分量qx,qy,qz分別為

其中,R3為三維實數(shù)空間,dW=dVxdVydVz(下同),m為分子質(zhì)量,n為分子數(shù)密度,kB為玻爾茲曼常數(shù);i,j=x,y,z;C為分子熱運動速度,Ci=Vi?Ui,etr為單位質(zhì)量動量.此外流場壓力p,黏性切應(yīng)力張量分量 τij分別為

這里,符號 δij表示當i=j時為1,當i≠j時為0.

氣體動理論統(tǒng)一算法所用玻爾茲曼方程可計算模型[24-31]可以寫為

這里,ν 為氣體分子碰撞松弛參數(shù),fN為當?shù)仄胶鈶B(tài)速度分布函數(shù).常用的沙克霍夫模型[12]和貝利模型[15]的表達式分別為

這里fM為麥克斯韋速度分布函數(shù),即

任何一個模型方程都必須滿足玻爾茲曼方程碰撞積分的如下3 個基本屬性[39].

(1)守恒性(質(zhì)量、動量、能量守恒),即

(3)在平衡態(tài)時,有Sm(fM)=0.

此外,對于單原子氣體,模型方程的查普曼?恩斯科分析得到的普朗特數(shù)Pr應(yīng)約為2/3.

為研究構(gòu)造新的玻爾茲曼方程可計算模型,首先分析其一階查普曼?恩斯科分析解[9],其可以寫為

顯然,Φ (C)的第一項與常用的沙克霍夫模型[12]中第二項除系數(shù)外表達式一致,而 Φ (C)的第二項與貝利模型[15]的第二項除系數(shù)外表達式也一致.

從玻爾茲曼方程的查普曼?恩斯科分析解可以看出,流場中氣體的速度分布函數(shù)與局部宏觀流動參數(shù)如速度、熱流、應(yīng)力等相關(guān),而沙克霍夫模型沒有考慮應(yīng)力的影響,貝利模型沒有考慮熱流的影響.當流場中局部區(qū)域的溫度、速度梯度較大即熱流和應(yīng)力影響顯著時,將出現(xiàn)較為嚴重的非平衡效應(yīng),僅考慮單一因素難以準確模擬復(fù)雜非平衡現(xiàn)象.因此,可以依據(jù)玻爾茲曼方程的查普曼?恩斯科分析解以及常用模型方程的表達式,構(gòu)造如下同時考慮熱流和應(yīng)力影響的新型玻爾茲曼方程可計算模型

這里,A和B是與分子速度無關(guān)的常量,νN為碰撞松弛時間的倒數(shù).顯然,當A=1?Pr,B=0時對應(yīng)沙克霍夫模型[12],A=0,B=(Pr?1)/Pr時對應(yīng)貝利模 型[15].

2 新型玻爾茲曼模型方程的數(shù)學(xué)分析

在構(gòu)造玻爾茲曼方程的一個新模型后,需要對其性質(zhì)進行分析,驗證其是否滿足玻爾茲曼方程碰撞積分項的基本屬性.

2.1 守恒性

由于

其中h=etr+p/ρ.上述3 個方程分別對應(yīng)質(zhì)量、動量和能量守恒方程.由此可知新模型方程滿足玻爾茲曼方程碰撞積分項的守恒性.

在平衡態(tài)時熱流和應(yīng)力張量均為零,顯然此時當?shù)佧溈怂鬼f分布是模型方程的解.

2.2 查普曼?恩斯科分析

為開展查普曼?恩斯科分析,首先將新模型方程的無量綱形式寫成如下形式

其中 ξ 為與來流克努森數(shù)(Kn)同量級的小量,頂標 ∧表示無量綱量.按查普曼?恩斯科分析的流程[9],定義

進而黏性系數(shù) μ、熱傳導(dǎo)系數(shù) κ、普朗特數(shù)Pr為

可知對于沙克霍夫模型[12],A=1?Pr,B=0,νS=p/μ ;對于貝利模型[15],A=0,B=(Pr?1)/Pr,νVVB=Pr·p/μ,這與前述一致.通過選擇合適的A,B可以得到正確的普朗特數(shù)Pr.

2.3 高階矩分析

為了進一步驗證新模型方程在數(shù)學(xué)上與玻爾茲曼方程的相容性,這里從麥克斯韋分子滿足的玻爾茲曼方程高階矩進行分析.由于麥克斯韋分子的特殊性,即其碰撞積分與相對速度無關(guān),可以得到對應(yīng)的高階矩.

定義碰撞積分項的矩為

該結(jié)論與上一節(jié)一致,這進一步證明了新模型方程與玻爾茲曼方程的相容性,與沙克霍夫模型[12]和貝利模型[15]的一致性.

2.4 H 定理的證明

定義[39]

且從查普曼?恩斯科分析的過程可知分布函數(shù)f相對fM的擾動為一小量,因此有

即當上述不等式成立,且Kn較小時,新模型方程滿足H 定理.從上面數(shù)學(xué)分析過程來看,νcoe在其取值范圍內(nèi)可以任意選取,但在實際物理過程中即分布函 數(shù)的碰撞松弛中 νcoe需要具有一定的物理意義.

2.5 碰撞松弛參數(shù)的確定

在氣體動理學(xué)理論中,根據(jù)碰撞間隔理論[39]可以推導(dǎo)出碰撞時間間隔 η、黏性系數(shù)μ、壓力p之間的關(guān)系為

顯然,若將 η 直接視為玻爾茲曼模型方程中氣體分子碰撞松弛參數(shù) ν 的倒數(shù),則上式與沙克霍夫模型一致,對應(yīng)于 νcoe=1.但上述碰撞間隔沒有考慮到分子碰撞后的速度殘留[39],因此需要對其進行修正,通常將其取為擾動馳豫時間,即

表征為分布函數(shù)松弛至其初始值的 1/e所歷經(jīng)的時間(e 為自然對數(shù)的底),這里 λ 為平均自由程.

此外在直接模擬蒙特卡洛方法中,對于變徑硬球(VHS)模型有[41]

上述 νcoe是通過微觀分子碰撞間隔理論得到的,在氣體動理學(xué)介觀理論中,為了表征分布函數(shù)f松弛至的1/e 的物理意義,同時使 νcoe滿足式(60),本文針對新模型方程取

對于單原子氣體氬氣(Ar),ω=0.81,此時νcoe=0 .819;對于空氣,ω=0.77,νcoe=0.788 4.

3 數(shù)值分析

3.1 一維激波結(jié)構(gòu)模擬分析

為驗證本文新建模型的有效性,首先模擬了一維激波結(jié)構(gòu).設(shè)置激波馬赫數(shù)為Mas=8,上游條件為:數(shù)密度n1=1 × 1023m?3,溫度T1=200 K,速度V1=2105.738 m/s;下游條件:n2=3.821 × 1020m?3,T2=4174.414 K,V2=551.111 m/s;氣體為氬氣;計算中激波上下游取當?shù)仄胶夥植?空間網(wǎng)格間距取上游氣體分子平均自由程的0.05 倍.

圖1 給出了采用沙克霍夫模型、本文新建模型、貝利模型模擬的激波輪廓線與直接模擬蒙特卡洛方法結(jié)果的對比,分別用“Shakhov”“new model”“VVB”和“DSMC”表示,其中直接模擬蒙特卡洛方法結(jié)果由貝德(Bird)的一維可視化直接模擬(DS1V)代碼[41]計算得到,λ1表示激波上游氣體子平均自由程,并將密度輪廓 (ρ?ρ1)/(ρ2?ρ1)=0.5 置于x/λ1=0 處.圖2 給出了3 種模型模擬的一維激波內(nèi)部無量綱熱流和切應(yīng)力變化曲線.從圖1 可以看出本文新建模型能較好地反映Mas=8 的激波輪廓,驗證了其有效性;對于幾種模型而言,沙克霍夫模型與直接模擬蒙特卡洛方法結(jié)果符合最好,本文模型次之;差別主要在激波上游靠近x/λ1=0 區(qū)域,從圖2 可以看出該區(qū)域無量綱切應(yīng)力 τxx變化很小,而無量綱熱流qx變化較大、且各模型結(jié)果區(qū)別也較大,僅考慮了熱流項的沙克霍夫模型與直接模擬蒙特卡洛方法結(jié)果最為接近,說明在一維激波內(nèi)部宏觀切應(yīng)力對內(nèi)部參數(shù)分布影響很小,這可能是由于一維問題的宏觀切應(yīng)力僅出現(xiàn)在x方向且無量綱值相對熱流較小,黏性效應(yīng)不顯著.此外模型方程模擬的溫度輪廓在激波上游均出現(xiàn)了提前抬起的現(xiàn)象,這 可能是碰撞頻率與分子速度無關(guān)造成的.

圖1 不同模型模擬的激波輪廓與直接模擬蒙特卡洛方法對比Fig.1 Comparison of shock profiles simulated by GKUA with different models and DSMC

圖2 不同模型模擬的無量綱熱流和切應(yīng)力對比Fig.2 Comparison of non-dimensional heat fluxes and shear stresses simulated by GKUA with different models

3.2 平板繞流模擬分析

為研究新模型方程模擬飛行器再入時在近空間環(huán)境黏性效應(yīng)急劇增加的非平衡流場中黏性應(yīng)力和熱流影響較大的區(qū)域流動問題,首先考察平板繞流.通常對于這類問題,在克努森層內(nèi)越靠近平板前緣,溫度和速度梯度越大,非平衡效應(yīng)越顯著,對應(yīng)的熱流和黏性效應(yīng)影響也越發(fā)嚴重.

來流條件:馬赫數(shù)Ma∞=5,來流氣體為氬氣,來流溫度和物面溫度均為200 K,數(shù)密度為n∞=1.175 × 1020m?3.平板長度為1 m、寬0.04 m.計算中,為了揭示高馬赫數(shù)附面層內(nèi)流動與傳熱機制,直接模擬蒙特卡洛方法的背景網(wǎng)格間距取3.62 mm,統(tǒng)一算法貼近物面第一排網(wǎng)格間距為2.5 mm,均遠小于來流分子平均自由程10 mm.統(tǒng)一算法模擬中除了所用新、舊碰撞模型外其他計算條件完全一樣,物面采用與直接模擬蒙特卡洛方法一致的無穿透麥克斯韋散射分布,來流邊界為來流無量綱麥克斯韋平衡分布,出口邊界根據(jù)當?shù)胤肿舆\動方向按特征邊界條件處理.圖3 給出了本文采用新模型方程模擬的溫度分布云圖(用GKUA_NewModel 表示的下半流場)與二維可視化直接模擬(DS2V)軟件[42]模擬結(jié)果(用DSMC 表示的上半流場)的對比情況.從圖中可以看出,兩個結(jié)果幾乎完全一致,而直接模擬蒙特卡洛方法尤其在尾部低速流場存在明顯的統(tǒng)計波動.這說明新模型模擬的流場結(jié)果是可信的.

圖3 本文模擬的溫度云圖與直接模擬蒙特卡洛方法結(jié)果對比Fig.3 Comparison of temperature contours simulated by GKUA with new model and DSMC

圖4 給出了在x=0.1 m 截面上壓力和溫度沿y方向不同碰撞模型的統(tǒng)一算法模擬結(jié)果與直接模擬蒙特卡洛方法的對比,其中本文建立的新模型用“new model”表示,沙克霍夫模型[12]用“Shakhov”表示,貝利模型[15]用“VVB”表示.從圖中可以看出,整體上使用新模型的統(tǒng)一算法結(jié)果與直接模擬蒙特卡洛方法模擬值相容,局部附面層區(qū)域y<0.15 m 時,相比沙克霍夫模型,新模型統(tǒng)一算法結(jié)果與直接模擬蒙特卡洛方法模擬值較為接近;而且還可看出對此黏性應(yīng)力影響極為嚴重的附面層流動,使用貝利模型的統(tǒng)一算法結(jié)果與直接模擬蒙特卡洛方法模擬值符合更好;在0.15 m 0.225 m)流動處于近平衡態(tài),直接模擬蒙特卡洛方法依靠麥克斯韋平衡態(tài)分布隨機采樣模擬分子的運動與碰撞,而統(tǒng)一算法基于求解當?shù)貧怏w分子速度分布函數(shù)的演化更新實時捕捉宏觀流動量變化分布,使得不同模型得到的統(tǒng)一算法計算值較為一致,均與直接模擬蒙特卡洛方法模擬存在一定的偏差,反映出進一步發(fā)展求解玻爾茲曼模型速度分布函數(shù)方程的氣體動理論統(tǒng)一算法有其獨到新穎性與應(yīng)用價值.

圖4 x=0.1 m 截面上壓力和溫度沿y 方向不同碰撞模型的統(tǒng)一算法計算值與直接模擬蒙特卡洛方法模擬結(jié)果的對比Fig.4 Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of x=0.1 m

圖5 和圖6 分別給出了在x=0.5 m 和0.8 m 截面上壓力和溫度沿y方向不同碰撞模型的統(tǒng)一算法計算結(jié)果與直接模擬蒙特卡洛方法模擬值對比.隨著x增大剪切層范圍也擴大,在這些區(qū)域內(nèi)流體黏性的影響較大,因而考慮了應(yīng)力張量的碰撞模型,如使用3 種碰撞模型得到的統(tǒng)一算法計算分布變化趨勢彼此一致,尤其是貝利模型與本文新提出的模型方程統(tǒng)一算法計算分布近乎完全重合,能較好模擬這些區(qū)域的局部流動,整體上都與直接模擬蒙特卡洛方法模擬結(jié)果吻合.而當y增大到剪切層外緣時,熱流的影響開始顯現(xiàn)加劇的某些區(qū)域,使用沙克霍夫模型得到的統(tǒng)一算法結(jié)果更吻合直接模擬蒙特卡洛方法模擬值.同樣在當y進一步增大到遠離附面層臨近計算區(qū)域外緣,流動梯度仍然較大的近平衡態(tài)局部流場,幾種模型得到的統(tǒng)一算法計算值均與直接模擬蒙特卡洛方法模擬值存在一定的系統(tǒng)誤差,究其實質(zhì)極可能緣于直接模擬蒙特卡洛方法完全使用麥克斯韋平衡態(tài)分布隨機統(tǒng)計模擬且存在統(tǒng)計漲落,導(dǎo)致了兩種方法體制的完全不一樣,對于流場梯度較大的局部流動模擬,彼此結(jié)果間存在一定的計算偏差.

圖5 x=0.5 m 截面上壓力和溫度沿y 方向不同碰撞模型統(tǒng)一算法計算結(jié)果與直接模擬蒙特卡洛方法模擬值比較Fig.5 Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of x=0.5 m

圖6 x=0.8 m 截面上壓力和溫度沿y 方向不同碰撞模型統(tǒng)一算法模擬結(jié)果與直接模擬蒙特卡洛方法的對比Fig.6 Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of x=0.8 m

從圖4～ 圖6 可以看出,由于本文建立的玻爾茲曼方程新模型同時兼顧了流場中應(yīng)力和熱流對氣體分子速度分布函數(shù)的影響,整體而言能較好地模擬全域流動.從統(tǒng)一算法計算x=0.1 m,0.5 m 和0.8 m 壓力和溫度變化來看,物面壓力從1.5 Pa 降至0.9 Pa和0.8 Pa,溫度從380 K 降至320 K 和300 K,相對于壁溫200 K,物面溫度跳躍逐漸降低,說明高馬赫數(shù)繞流平板前緣強壓縮流動在近空間飛行環(huán)境稀薄效應(yīng)影響下逐漸變得稀疏,沿平板物面壓力、溫度不斷減小,過渡到平板尾部低壓常規(guī)流.

圖7 給出了駐點線上不同碰撞模型得到的統(tǒng)一算法計算溫度分布與直接模擬蒙特卡洛方法模擬值對比.平板頭部壓縮激波層是非平衡效應(yīng)最為嚴重的區(qū)域,溫度梯度和速度梯度均較大.從圖中可以看出,在近空間飛行環(huán)境激波內(nèi)部,本文建立的玻爾茲曼方程新模型統(tǒng)一算法計算結(jié)果始終與直接模擬蒙特卡洛方法模擬值符合較好.這進一步說明了需要同時考慮流場中黏性應(yīng)力和熱流對氣體分子速度分布函數(shù)演化更新的影響.同時也看出基于玻爾茲曼方程可計算建模統(tǒng)一算法計算結(jié)果在駐點激波前部流動趨于近平衡態(tài)時溫度高于直接模擬蒙特卡洛方法模擬值,這與一維激波結(jié)構(gòu)中溫度輪廓的變化規(guī)律類似.

圖7 駐點線上不同碰撞模型統(tǒng)一算法計算的溫度分布與直接模擬蒙特卡洛方法模擬值比較Fig.7 Comparison of temperature simulated by GKUA with different collision models and DSMC along the stagnation line

3.3 雙圓柱干擾流動模擬分析

為驗證新建模型統(tǒng)一算法對激波干擾的模擬能力,這里考察雙圓柱干擾流動.來流條件和計算條件及邊界條件與上節(jié)一致,圓柱直徑為1 m,兩圓柱圓心間距為4 m.在該條件下上下兩個圓柱的頭部脫體激波會相互作用,形成復(fù)雜的激波干擾流動.

圖8 給出了分別使用3 種模型方程的統(tǒng)一算法計算溫度分布云圖與直接模擬蒙特卡洛方法模擬值的對比.兩圓柱脫體激波相互作用后在x=?0.19m 的位置形成一個正激波類馬赫桿,跨過馬赫桿后溫度和壓力增加、速度降低.而馬赫桿兩端形成的斜激波僅對圓柱尾流產(chǎn)生影響,將導(dǎo)致單邊圓柱法向力很小.從云圖可以看出,使用貝利模型統(tǒng)一算法得到核心區(qū)溫度分布與直接模擬蒙特卡洛方法模擬值最為接近,但馬赫桿位置稍向后,而沙克霍夫模型差別最大,對馬赫桿位置而言本文新建模型統(tǒng)一算法與直接模擬蒙特卡洛方法符合最好.表1 給出了3 種碰撞模型統(tǒng)一算法及直接模擬蒙特卡洛方法模擬得到的單邊圓柱軸向力系數(shù)及法向力系數(shù),以及不同模型所得軸向力系數(shù)與直接模擬蒙特卡洛方法模擬結(jié)果的偏差.從表中可以看出,本文新建立的碰撞模型統(tǒng)一算法計算結(jié)果與直接模擬蒙特卡洛方法模擬值最大偏差在4.6%以內(nèi),能較好地模擬氣動力特性,精度優(yōu)于常用的沙克霍夫模型.

圖8 3 種模型方程計算得到的溫度分布云圖與直接模擬蒙特卡洛方法結(jié)果的對比Fig.8 Comparison of temperature contours simulated by GKUA with three collision models and DSMC

表1 不同模型/方法所得軸向力系數(shù)及其偏差與法向力系數(shù)Table 1 The axial and normal force coefficients by different models and DSMC method and errors between axial force coefficients

圖9 給出了對稱面(即y=0)上馬赫桿附近沿x方向不同碰撞模型統(tǒng)一算法計算得到的壓力、溫度和速度分布與直接模擬蒙特卡洛方法模擬值的對比情況.本文新建碰撞模型較好地求解了二維繞流的激波位置及激波前后宏觀參數(shù)的變化趨勢,明顯優(yōu)于沙克霍夫模型,在激波內(nèi)部模擬結(jié)果優(yōu)于貝利模型.這是因為激波前后溫度和速度梯度較大,但y=0 截面上僅存在x方向熱流qx,而x和y方向的切應(yīng)力 τxx和 τyy同時存在,速度梯度對分布函數(shù)的影響更大,使得本文新建碰撞模型能更好地模擬激波位置和內(nèi)部參數(shù)變化.這與模擬一維激波結(jié)構(gòu)所得結(jié)論不一致,可能是由于隨著維度的增加多方向宏觀切應(yīng)力導(dǎo)致黏性效應(yīng)更加顯著,其對分子速度分布函數(shù)的影響增大,與熱流共同作用影響激波位置 及其內(nèi)部參數(shù)的分布.

圖9 y=0 m 截面上壓力、溫度和速度沿 x 方向不同碰撞模型模擬結(jié)果與直接模擬蒙特卡洛方法的對比Fig.9 Comparison of pressure,temperature and velocity simulated by GKUA with different collision models and DSMC along the direction of x on the section of y=0 m

4 結(jié)論

本文通過分析玻爾茲曼方程的一階查普曼?恩斯科分析解,結(jié)合沙克霍夫模型和貝利模型構(gòu)造過程及數(shù)學(xué)推導(dǎo)與物理分析,建立了一種同時考慮熱流矢量和應(yīng)力張量影響、滿足玻爾茲曼方程高階碰撞矩的綜合可計算碰撞模型,并在數(shù)學(xué)上分析了其守恒性、H 定理等基本屬性,證明了新模型方程與玻爾茲曼方程的相容性,給出了新模型與現(xiàn)有模型的關(guān)系,通過碰撞動力學(xué)確定了碰撞松弛參數(shù)統(tǒng)一表達式.經(jīng)過模擬一維激波結(jié)構(gòu)、二維近空間飛行環(huán)境平板和多體圓柱干擾流動,并與直接模擬蒙特卡洛方法及不同碰撞模型方程的統(tǒng)一算法結(jié)果對比分析,得出以下結(jié)論:

(1)在稀薄效應(yīng)顯著、非平衡現(xiàn)象突出的流動區(qū)域,熱流矢量和應(yīng)力張量對當?shù)貧怏w分子速度分布函數(shù)的影響較大,在構(gòu)造碰撞松弛模型時需要同時考慮兩者的影響;

(2)在激波捕捉方面,流場中黏性效應(yīng)顯著的區(qū)域沙克霍夫模型結(jié)果與直接模擬蒙特卡洛方法模擬值差別較大,本文新建模型與貝利模型結(jié)果接近,在激波內(nèi)部本文模型與直接模擬蒙特卡洛方法符合更好,而在黏性效應(yīng)較小時僅考慮熱流時結(jié)果與直接模擬蒙特卡洛方法符合更好,進一步說需要考慮多參數(shù)對分子速度分布函數(shù)的綜合影響;

(3)本文新建模型的碰撞松弛參數(shù)中引入了與自然對數(shù)的底 e 相關(guān)的常數(shù)作為自由參數(shù),間接表征了分布函數(shù)松弛的物理意義,雖然能滿足模型方程的H 定理,但在流動的某些區(qū)域仍然與直接模擬蒙特卡洛方法存在差距.若能構(gòu)造一個與當?shù)胤瞧胶鈪?shù)相關(guān)聯(lián)的碰撞松弛參數(shù),將氣體分子速度分布函數(shù)的松弛過程與當?shù)亓鲃訁?shù)的梯度耦合,應(yīng)該能更好地模擬再入跨流域非平衡流動過程,這也是本文下一步擬開展的工作.