一種新的玻爾茲曼方程可計(jì)算模型構(gòu)造與分析1)

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

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

? (北京航空航天大學(xué)國家計(jì)算流體力學(xué)實(shí)驗(yàn)室,北京 100191)

引言

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

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

模型方程作為玻爾茲曼方程的近似,必須滿足玻爾茲曼方程的基本屬性,如總和不變量性質(zhì)、H 定理、正確的輸運(yùn)系數(shù)等.第一個(gè)也是最著名的模型方程是巴特納?格羅斯?克魯克模型[8](BGK模型),但通過對其進(jìn)行查普曼?恩斯科(Chapman?Enskog)展開分析得到的普朗特?cái)?shù)(Pr)為1,而單原子氣體的Pr約為2/3[9].為克服這一缺陷,國內(nèi)外學(xué)者提出了一系列模型方程,包括橢球統(tǒng)計(jì)?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],上述模型方程均能滿足玻爾茲曼方程的基本屬性,且與原始玻爾茲曼方程的碰撞項(xiàng)具有類似的松弛特性,能保持與納維?斯托克斯方程、直接模擬蒙特卡洛方法的收斂一致,因結(jié)構(gòu)簡單、貼近跨流域流體物理變化特點(diǎn)、可計(jì)算技術(shù)針對性強(qiáng)而得到廣泛應(yīng)用[21],如泰塔雷夫(Titarev)建立了物理空間和速度空間均采用非結(jié)構(gòu)網(wǎng)格的內(nèi)斯韋塔依(Nesvetay)求解器[22-23],針對返回式空間飛行器等高馬赫數(shù)流動(dòng)與直接模擬蒙特卡洛方法進(jìn)行了系列對比分析.李志輝等[24-31]從研究描述各流域均適用的氣體分子速度分布函數(shù)方程出發(fā),吸收計(jì)算數(shù)學(xué)指數(shù)型積分求解原理,發(fā)展了氣體分子運(yùn)動(dòng)論離散速度坐標(biāo)法,研制經(jīng)改進(jìn)的高斯?埃爾米特積分方法等系列離散速度數(shù)值積分技術(shù),消除氣體分子速度分布函數(shù)對速度空間的連續(xù)依賴性,將玻爾茲曼模型方程化為在各個(gè)離散速度坐標(biāo)點(diǎn)處具有非線性源項(xiàng)的雙曲型守恒方程,基于計(jì)算流體力學(xué)數(shù)值求解技術(shù),構(gòu)造直接求解速度分布函數(shù)演化更新的氣體動(dòng)理論數(shù)值格式和大規(guī)模并行計(jì)算策略,建立了氣體動(dòng)理論統(tǒng)一算法(GKUA),開展了自由分子流到連續(xù)流復(fù)雜高超聲速氣動(dòng)力/熱問題大規(guī)模并行計(jì)算及其在航天工程中的應(yīng)用研究,同時(shí)對跨流域多體干擾流動(dòng)問題、跨流域非定常流動(dòng)問題進(jìn)行了研究,獲得較好的成功.徐昆等[32-37]利用上述離散速度坐標(biāo)法與BGK 類型格式相耦合,發(fā)展了統(tǒng)一氣體動(dòng)理學(xué)格式(UGKS),使之適于模擬稀薄流到連續(xù)流氣體流動(dòng)問題.

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

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

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

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

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

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

這里,符號(hào) δij表示當(dāng)i=j時(shí)為1,當(dāng)i≠j時(shí)為0.

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 守恒性

由于

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

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

2.2 查普曼?恩斯科分析

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

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

進(jìn)而黏性系數(shù) μ、熱傳導(dǎo)系數(shù) κ、普朗特?cái)?shù)Pr為

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

2.3 高階矩分析

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

定義碰撞積分項(xiàng)的矩為

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

2.4 H 定理的證明

定義[39]

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

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

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

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

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

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

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

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

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

3 數(shù)值分析

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

為驗(yàn)證本文新建模型的有效性,首先模擬了一維激波結(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;氣體為氬氣;計(jì)算中激波上下游取當(dāng)?shù)仄胶夥植?空間網(wǎng)格間距取上游氣體分子平均自由程的0.05 倍.

圖1 給出了采用沙克霍夫模型、本文新建模型、貝利模型模擬的激波輪廓線與直接模擬蒙特卡洛方法結(jié)果的對比,分別用“Shakhov”“new model”“VVB”和“DSMC”表示,其中直接模擬蒙特卡洛方法結(jié)果由貝德(Bird)的一維可視化直接模擬(DS1V)代碼[41]計(jì)算得到,λ1表示激波上游氣體子平均自由程,并將密度輪廓 (ρ?ρ1)/(ρ2?ρ1)=0.5 置于x/λ1=0 處.圖2 給出了3 種模型模擬的一維激波內(nèi)部無量綱熱流和切應(yīng)力變化曲線.從圖1 可以看出本文新建模型能較好地反映Mas=8 的激波輪廓,驗(yàn)證了其有效性;對于幾種模型而言,沙克霍夫模型與直接模擬蒙特卡洛方法結(jié)果符合最好,本文模型次之;差別主要在激波上游靠近x/λ1=0 區(qū)域,從圖2 可以看出該區(qū)域無量綱切應(yīng)力 τxx變化很小,而無量綱熱流qx變化較大、且各模型結(jié)果區(qū)別也較大,僅考慮了熱流項(xiàng)的沙克霍夫模型與直接模擬蒙特卡洛方法結(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 平板繞流模擬分析

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

來流條件:馬赫數(shù)Ma∞=5,來流氣體為氬氣,來流溫度和物面溫度均為200 K,數(shù)密度為n∞=1.175 × 1020m?3.平板長度為1 m、寬0.04 m.計(jì)算中,為了揭示高馬赫數(shù)附面層內(nèi)流動(dòng)與傳熱機(jī)制,直接模擬蒙特卡洛方法的背景網(wǎng)格間距取3.62 mm,統(tǒng)一算法貼近物面第一排網(wǎng)格間距為2.5 mm,均遠(yuǎn)小于來流分子平均自由程10 mm.統(tǒng)一算法模擬中除了所用新、舊碰撞模型外其他計(jì)算條件完全一樣,物面采用與直接模擬蒙特卡洛方法一致的無穿透麥克斯韋散射分布,來流邊界為來流無量綱麥克斯韋平衡分布,出口邊界根據(jù)當(dāng)?shù)胤肿舆\(yùn)動(dòng)方向按特征邊界條件處理.圖3 給出了本文采用新模型方程模擬的溫度分布云圖(用GKUA_NewModel 表示的下半流場)與二維可視化直接模擬(DS2V)軟件[42]模擬結(jié)果(用DSMC 表示的上半流場)的對比情況.從圖中可以看出,兩個(gè)結(jié)果幾乎完全一致,而直接模擬蒙特卡洛方法尤其在尾部低速流場存在明顯的統(tǒng)計(jì)波動(dò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 時(shí),相比沙克霍夫模型,新模型統(tǒng)一算法結(jié)果與直接模擬蒙特卡洛方法模擬值較為接近;而且還可看出對此黏性應(yīng)力影響極為嚴(yán)重的附面層流動(dòng),使用貝利模型的統(tǒng)一算法結(jié)果與直接模擬蒙特卡洛方法模擬值符合更好;在0.15 m 0.225 m)流動(dòng)處于近平衡態(tài),直接模擬蒙特卡洛方法依靠麥克斯韋平衡態(tài)分布隨機(jī)采樣模擬分子的運(yùn)動(dòng)與碰撞,而統(tǒng)一算法基于求解當(dāng)?shù)貧怏w分子速度分布函數(shù)的演化更新實(shí)時(shí)捕捉宏觀流動(dòng)量變化分布,使得不同模型得到的統(tǒng)一算法計(jì)算值較為一致,均與直接模擬蒙特卡洛方法模擬存在一定的偏差,反映出進(jìn)一步發(fā)展求解玻爾茲曼模型速度分布函數(shù)方程的氣體動(dòng)理論統(tǒng)一算法有其獨(dú)到新穎性與應(yīng)用價(jià)值.

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

圖5 x=0.5 m 截面上壓力和溫度沿y 方向不同碰撞模型統(tǒng)一算法計(jì)算結(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 可以看出,由于本文建立的玻爾茲曼方程新模型同時(shí)兼顧了流場中應(yīng)力和熱流對氣體分子速度分布函數(shù)的影響,整體而言能較好地模擬全域流動(dòng).從統(tǒng)一算法計(jì)算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ù)繞流平板前緣強(qiáng)壓縮流動(dòng)在近空間飛行環(huán)境稀薄效應(yīng)影響下逐漸變得稀疏,沿平板物面壓力、溫度不斷減小,過渡到平板尾部低壓常規(guī)流.

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

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

3.3 雙圓柱干擾流動(dòng)模擬分析

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

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

圖8 3 種模型方程計(jì)算得到的溫度分布云圖與直接模擬蒙特卡洛方法結(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)一算法計(jì)算得到的壓力、溫度和速度分布與直接模擬蒙特卡洛方法模擬值的對比情況.本文新建碰撞模型較好地求解了二維繞流的激波位置及激波前后宏觀參數(shù)的變化趨勢,明顯優(yōu)于沙克霍夫模型,在激波內(nèi)部模擬結(jié)果優(yōu)于貝利模型.這是因?yàn)榧げㄇ昂鬁囟群退俣忍荻容^大,但y=0 截面上僅存在x方向熱流qx,而x和y方向的切應(yīng)力 τxx和 τyy同時(shí)存在,速度梯度對分布函數(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)與物理分析,建立了一種同時(shí)考慮熱流矢量和應(yīng)力張量影響、滿足玻爾茲曼方程高階碰撞矩的綜合可計(jì)算碰撞模型,并在數(shù)學(xué)上分析了其守恒性、H 定理等基本屬性,證明了新模型方程與玻爾茲曼方程的相容性,給出了新模型與現(xiàn)有模型的關(guān)系,通過碰撞動(dòng)力學(xué)確定了碰撞松弛參數(shù)統(tǒng)一表達(dá)式.經(jīng)過模擬一維激波結(jié)構(gòu)、二維近空間飛行環(huán)境平板和多體圓柱干擾流動(dòng),并與直接模擬蒙特卡洛方法及不同碰撞模型方程的統(tǒng)一算法結(jié)果對比分析,得出以下結(jié)論:

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

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

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