皮興才，李志輝*，彭傲平，張子彬

（1.中國空氣動力研究與發(fā)展中心超高速空氣動力研究所，綿陽621000；2.國家計算流體力學實驗室，北京100191）

1 引言

服役期滿的航天器將再入、解體并進入近空間飛行環(huán)境。由于稠密大氣對解體物的氣動力、熱作用，大部解體物會在再入過程中熔融、燒蝕、蒸發(fā)消亡掉。但隨著飛行高度、速度、能量不斷降低，至少10%～40%質(zhì)量的殘骸碎片存留，形成板、方/圓柱、球形物等不同類型的殘骸，隕落到地面，可能對人財物、生態(tài)系統(tǒng)安全造成威脅［1-2］。

格子Boltzmann方法（Lattice Boltzmannn Method，LBM）起源于格子氣自動機方法，由Mc-Namara和Zanetti提出［3］，并由Higuera［4］將其應用到實際的流體力學數(shù)值計算中。目前已經(jīng)成為一種新的計算流體力學及數(shù)值傳熱學分析方法，被廣泛應用于多孔介質(zhì)流動、多相流動、微流動等問題分析。

為了保證離散守恒性，標準LBM的平衡態(tài)分布函數(shù)由Boltzmann方程的經(jīng)典解析解：Maxwell平衡態(tài)速度分布函數(shù)作小Ma數(shù)Taylor展開獲得。這使得通過Chapman-Enskog展開標準LBM獲得的宏觀流動物理量方程是具有Ma2量級的附加項的不可壓縮Navier-Stokes方程［5］。因此，標準LBM被視為求解不可壓縮流動的壓力修正方法［6］。為了使格子Boltzmann方法能分析可壓縮流動，Shi等［7］提出的求解可壓縮Euler方程的格子Boltzmann模型，Chen等［8］提出的求解可壓縮Navier-Stokes方程的格子Boltzmann模型。但這些模型無法實現(xiàn)比熱比γ的可調(diào)，且其還原的比熱比不滿足實際物理意義。Kataoka等［9］提出的多速度格子模型克服了上述問題，并被成功應用到Ma小于1的可壓縮粘性流動模擬。

為了使LBM方法具備高馬赫數(shù)可壓縮流動的模擬能力，Qu等［10］提出用圓函數(shù)替代經(jīng)典的Maxwell平衡態(tài)分布函數(shù)來構(gòu)造格子Boltzmann離散速度模型，并確保圓函數(shù)滿足恢復可壓縮N-S方程所需要的矩關(guān)系，成功應用于可壓縮流動模擬。這樣的模型構(gòu)造思路吸引了眾多研究者：如Li等［11］將其與耦合雙分布函數(shù)（Coupled Double Distribution Function，CDDF）Boltzmann模型方程結(jié)合，成功應用于Riemann問題、雙馬赫反射、Couette流等經(jīng)典流動問題的模擬，并實現(xiàn)了普朗特數(shù)Pr、比熱比γ任意可調(diào)；Qiu等［12］將該方法成功應用于可壓縮邊界層、激波邊界層干擾等流動模擬；Gan等［13-14］采用高階離散平衡態(tài)分布函數(shù)模型，構(gòu)造了滿足更復雜動理學矩關(guān)系的離散Boltzmann模型，成功將格子Boltzmann方法由可壓縮連續(xù)介質(zhì)流域模擬擴展到Burnett層面的可壓縮近連續(xù)流模擬，并揭示了更具普適性本構(gòu)關(guān)系，為高速可壓縮流動中的熱力學非平衡效應研究提供了新的建模思路。

綜上可看出：LBM在求解可壓縮流動方面，已從發(fā)展之初的格子氣思想，逐漸過渡到利用Boltzmann介觀思想構(gòu)造求解流體物理問題的數(shù)學關(guān)系，并采用格子速度模型進行模型構(gòu)建階段。但這些可壓縮格子模型在高馬赫數(shù)流動中應用案例還鮮有體現(xiàn)。為將LBM推廣應用于服役期滿航天器再入解體殘骸碎片的可壓縮流動問題，本文將構(gòu)造有限體積格子Boltzmann數(shù)值格式，求解耦合雙分布函數(shù)模型方程以表征非平衡屬性，并采用以圓函數(shù)為基礎發(fā)展的D2Q13離散速度模型表征平衡態(tài)，由此建立適于近空間連續(xù)流區(qū)非規(guī)則殘骸碎片繞流模擬的有限體積隱式格子Boltzmann方法，并在對Riemann問題、平板雙馬赫反射、RAE2822翼型跨聲速繞流模擬驗證基礎上，開展近空間多次解體形成的殘骸碎片簡化方柱模型超聲速繞流計算及其與N-S方程解算器結(jié)果對比檢驗，分析經(jīng)改進的格子Boltzmann理論模型在模擬近空間連續(xù)流區(qū)殘骸碎片繞流模擬潛力。

2 理論模型

因耦合雙分布函數(shù)Boltzmann模型方程［11，15］的普朗特數(shù)Pr、比熱比γ任意可調(diào)，且具有良好的數(shù)值穩(wěn)定性，本文采用它作為控制方程。其密度及能量分布函數(shù)的演化方程如式（1）所示：

式中，下標α代表不同的離散速度編號，fα是密度分布函數(shù)，hα是能量分布函數(shù)。是平衡態(tài)密度分布函數(shù)，是平衡態(tài)能量分布函數(shù)。eα是離散速度矢量，u是宏觀速度矢量，τf與τh分別是密度分布函數(shù)松弛時間及能量分布松弛時間，其滿足關(guān)系式（2）：

普朗特數(shù)Pr與松弛時間的關(guān)系為式（2）：

為了使演化方程通過Chapman-Enskog多尺度展開后能夠還原到可壓縮N-S方程，平衡態(tài)密度分布函數(shù)及平衡態(tài)能量分布函數(shù)需要滿足如式（4）所示關(guān)系式：

為了滿足上述關(guān)系，本文采用基于Qu提出的圓函數(shù)重構(gòu)平衡態(tài)分布函數(shù)的D2Q13離散速度模型［10-11］。D2Q13離散速度模型如式（5）所示：

其他離散速度點的平衡態(tài)密度分布函數(shù)滿足式（7）：

各離散速度點的平衡態(tài)能量分布函數(shù)見式（8）：

3 數(shù)值方法

前期應用Boltzmann方程可計算建模氣體動理論統(tǒng)一算法，求解服役期滿航天器解體殘骸的跨流域繞流問題表明：隨著飛行高度降低，解體殘骸碎片再凝固形成的非規(guī)則解體物速度迅速降低，呈現(xiàn)近地面連續(xù)/近連續(xù)流區(qū)非規(guī)則物形激波/膨脹波多物理場復雜繞流現(xiàn)象［16-17］。因此，引入有限體積隱式格式求解耦合雙分布函數(shù)格子Boltzmann模型方程，以保證構(gòu)造的數(shù)值格式具有更好的通量守恒性以及更強的非規(guī)則解體物外形適應能力。采用有限體積法處理演化方程式（1）：

式中，

可得式（11）：

式中，Ωij為結(jié)構(gòu)網(wǎng)格控制體 （i，j） 的體積，及表示F，S在單元控制體內(nèi)的均值，n為控制體界面沿i，j增大方向的法向單位向量。本文采用NND格式來重構(gòu)界面通量。因此，式（11）等號左邊第二項可以改寫為式（12）：

式中，Δs為控制體Ωij的邊長。利用Steger-Warming流通矢量分裂，將控制體界面上的通量分解為正通量和負通量，其中正通量由計

（15）：

在時間項的處理方面，本文采用多步龍格庫塔法（IMEX-RK）［6，18］。該方法采用顯式方式處理對流項，采用隱式方式處理源項，具有很強的剛性問題處理能力，被成功應用于BGK模型方程的求解。分布函數(shù)從第n時間步演化到n+1時間步的公式見式（16）：

式中，上標 ＜n＞表示第n步IMEX-RK變量，l表示IMEX-RK的級數(shù)。兩個l×l矩陣a~nn，ann以及l(fā)維向量w~n，wn為IMEX-RK方法的系數(shù)矩陣及系數(shù)向量。文獻［19］對IMEX-RK方法的理論及各階系數(shù)矩陣作了詳細的研討。

從式（16）可看出，第n步分布函數(shù)的計算需要第n步碰撞項的結(jié)果，其是典型的隱式格式。從演化方程式（9）可知，碰撞項包含分布函數(shù)、平衡態(tài)分布函數(shù)及碰撞松弛時間。由此可知，若能獲得第n步的宏觀物理量，此方程可顯式計算第n步分布函數(shù)，而不用迭代。值得注意的是Boltzmann模型方程具有質(zhì)量、動量及能量守恒性，即碰撞項對碰撞守恒量的矩為0。利用該性質(zhì)可獲得第n步的宏觀物理量：針對式中的密度分布函數(shù)分量，在等式各項分別乘以碰撞守恒量φ=（1，eα）T， 并 對 所 有 離 散 速 度 點 求 和 可 得 式（17）［11］：

同理，對能量分布函數(shù)進行處理可得式（18）：

由上式可看出，第n步宏觀物理量可通過前幾步的分布函數(shù)顯式地計算得出。有了各步的宏觀量，其對應的平衡態(tài)分布函數(shù)及碰撞松弛時間即可得到。

4 算例驗證與分析

采用上述方法進行Riemann問題、雙馬赫反射等經(jīng)典問題的模擬，以驗證算法程序可靠性，然后將此算法應用于翼型、方柱模型等近空間多次解體殘骸碎片簡化模型的亞跨超聲速繞流計算。從雙分布格子Boltzmann方程的構(gòu)造理念可看出，特征溫度的選取并不影響其還原物理量的速度分布函數(shù)矩關(guān)系。但特征溫度的選取直接影響離散速度點的平衡態(tài)分布函數(shù)的插值穩(wěn)定性，進而影響整個數(shù)值格式的穩(wěn)定性。經(jīng)驗表明，特征溫度的取值一般略大于流場總溫。

4.1 一維Riemann問題

一維Riemann問題流場包含激波、接觸間斷、膨脹波3類典型的可壓縮流動基本元素，流動具有解析解，并且由于流動的一維性，可方便可靠地進行邊界處理，使得考核能夠排除邊界條件對結(jié)果的影響，被廣泛用于代碼考核。本數(shù)值案例中，定義參考密度ρ0=1.165 kg/m3，參考溫度T0=303 K，參考粘性系數(shù)：μ0=1.86×10-5kg/( m·s)，比熱比γ=1.4，普朗特數(shù)Pr=0.72，參考速度u0以及參考壓力p0見式（19）：

案例1：SOD問題：

其中，L0=2m，特征時間t0=L0/u0，設定特征溫度Tc=2.5T0，時間步長設為d t=3000·μ0/p0，網(wǎng)格采用Nx×Ny=800×5的均勻網(wǎng)格，且d x=d y=L0/800。在x向邊界處，直接指定宏觀量，并由此得出邊界的平衡分布函數(shù)作為邊界條件；在y向，采用周期性邊界條件。圖1展示了在t=0.1644t0的數(shù)值計算結(jié)果及解析解的對比情況?？煽闯?，除激波后溫度場有稍許偏差外，數(shù)值結(jié)果整體與解析解吻合良好。

圖1 SOD問題數(shù)值解驗證Fig.1 Verification of numerical solution of SOD

案例2：強激波問題

這是一個典型的Colella爆炸波算例，其描述的是激波管隔膜破裂后的流動現(xiàn)象。由初始條件可看出，隔膜兩邊的壓力比大于10萬，對于數(shù)值方法的穩(wěn)定性、空間分辨力提出了巨大挑戰(zhàn)。本算例的基本參考物理量與上例保持一致，特征溫度Tc=1000T0。圖2展示了t=0.012t0的計算結(jié)果與解析解的對比情況。整體而言，計算結(jié)果與解析解吻合良好。當然，圖中也可看出壓力計算結(jié)果在激波后出現(xiàn)振蕩，進一步驗證發(fā)現(xiàn)不同的時間及空間格式會產(chǎn)生不同的振幅及頻率，并且一定程度的加密網(wǎng)格可以改善振蕩。密度以及溫度場在激波后與解析解出現(xiàn)稍許的偏離，研究認為這些偏離與空間離散格式有關(guān)［11，20］。

4.2 平板雙馬赫反射

服役期滿航天器再入解體過程高溫度梯度致物面附近流動壓縮氣動加熱積累，進入近空間連續(xù)流區(qū)多次解體殘骸碎片間繞流干擾出現(xiàn)壓縮激波反射嚴重的馬赫桿現(xiàn)象，激波的馬赫反射復雜流動一直是間斷激波多物理場用于理論數(shù)值方法檢驗的典型問題。為此，計算馬赫數(shù)為10的高超聲速來流以30°的角度入射平面后的流場結(jié)構(gòu)。本算例與一維Riemann問題數(shù)值模擬中采用的參考量保持一致。通過本算例，考核本文數(shù)值方法對二維可壓縮流動適應性，其物理模型如圖3所示。

計算區(qū)域為3 m×1 m的二維矩形域。反射面位于區(qū)域底面，反射面采用鏡面反射邊界條件，其余邊界的宏觀物理量根據(jù)運動激波關(guān)系式給出，密度及能量分布函數(shù)在流動入口邊界給定為平衡態(tài)，在出口邊界采用外推格式。計算中，特征溫度Tc=85T0。圖4給出t=0.062t0時刻的密度等值線及壓力等值線。由圖看出，計算結(jié)果清晰反映了激波馬赫桿，以及第二次反射激波及馬赫桿滑移線。還捕獲到馬赫數(shù)及斜劈角較大的情況下出現(xiàn)的第一馬赫桿突出變形現(xiàn)象，馬赫桿在接近壁面時逐漸垂直壁面的流動特點和理論分析吻合，與文獻［21，22］進行對比分析，也可發(fā)現(xiàn)相互吻合較好。

圖2 強激波問題數(shù)值解Fig.2 Numerical solution of strong shock

4.3 RAE2822翼型跨聲速繞流問題

RAE2822翼型是典型的二維跨音速流動算例，被16個EUROVAL歐洲合作項目組和AGARD挑選為算法驗證確認流動問題［23］。計算狀態(tài)為Ma=0.734，α=2.79°，Re=6.5×106，采用的網(wǎng)格由NPARC官網(wǎng)提供（如圖5），壁面采用非平衡外推方式實現(xiàn)無滑移速度邊界，無窮遠來流邊界條件按Maxwell平衡態(tài)分布給定，出口邊界條件由內(nèi)流場外推，特征溫度Tc=2.5T0。

圖3 雙馬赫反射問題的物理模型Fig.3 Physicalmodel of the double-M ach-reflection

圖4 雙馬赫反射流場密度、壓力等值線Fig.4 Density and pressure contours of the double-M ach-reflection

圖5 翼型繞流計算網(wǎng)格Fig.5 Com putationalmesh for flows around airfoil

圖6繪出該翼型面繞流壓力分布?？梢钥闯觯硇捅巢康臍饬鞅患铀俸?，壓力逐步降低。進一步，加速形成的超聲速氣流經(jīng)過背部激波后減速，壓力陡然回升。由于圖7給出了翼型表面壓力系數(shù)分布及其與實驗結(jié)果［24］對比情況?？梢钥闯?，除了背部湍流轉(zhuǎn)捩區(qū)及激波附著區(qū)之外，其他位置計算所得壓力系數(shù)與實驗結(jié)果吻合較好。

圖6 翼型繞流壓力等值線Fig.6 Pressure contours of the flows around airfoil

圖7 表面的壓力系數(shù)Fig.7 Pressure coefficient along the airfoil surface

4.4 解體殘骸碎片類方柱超聲速繞流問題

為考核本文算法在解體殘骸碎片類方柱超聲速繞流問題模擬能力，擬定二維方塊在來流條件飛行高度H=50 km的大氣環(huán)境下，以馬赫數(shù)Ma=2.0、雷諾數(shù)Re=8×104飛行的方塊繞流問題。計算中，特征溫度Tc=4.5T0，無窮遠處來流條件以Maxwell平衡態(tài)速度分布給定，出口條件由內(nèi)流場外推，壁面邊界條件采用鏡面反射滑移邊界條件。圖8展示了本文所發(fā)展的耦合雙分布函數(shù)有限體積格子Boltzmann方法的計算結(jié)果與N-S方程求解結(jié)果對比情況（“CDDF”表示耦合雙分布函數(shù)有限體積格子Boltz-mann方法結(jié)果，“NS”表示NS方法計算結(jié)果）。由圖可看出，繞流場密度、溫度、馬赫數(shù)分布與N-S方程計算結(jié)果整體吻合較好，但在溫度場的尾部回流區(qū)，本文計算結(jié)果（圖8（b））與N-S結(jié)果存在差異。進一步，采用等溫非平衡外推無滑移邊界條件計算結(jié)果的尾部回流區(qū)溫度場與N-S方程求解器的計算結(jié)果吻合良好，由此說明邊界條件對數(shù)值計算結(jié)果的合理性有明顯影響。隨著馬赫數(shù)增加，氣體分子速度分布函數(shù)的非平衡屬性顯現(xiàn)明顯，耦合雙分布函數(shù)有限體積格子Boltzmann方法的數(shù)值穩(wěn)定性變差，且尾部回流區(qū)域稀薄間斷粒子“真空”稀薄效應嚴重，會對數(shù)值穩(wěn)定性產(chǎn)生明顯影響。文獻［25］研究表明，采用更復雜的離散速度模型及多松弛時間演化方程可成功實現(xiàn)馬赫數(shù)5的前半圓柱繞流模擬，但具有回流區(qū)的高馬赫繞流模擬還鮮有體現(xiàn)，這也反映出格子Boltzmann方法對涉及非平衡效應基于速度分布函數(shù)高階取矩宏觀量的溫度、熱流量模擬存在局限性。對于近空間多次解體再凝固板、方/圓柱等殘骸/碎片近地面連續(xù)流，因飛行隕落速度很低，基本不存在非平衡效應，為此，本文發(fā)展經(jīng)改進耦合雙分布函數(shù)的有限體積格子Boltzmann方法對這類近空間連續(xù)/近連續(xù)流區(qū)殘骸碎片繞流問題模擬具有較好的適應性，進一步的分析確認有待今后開展。

5 結(jié)論

本文將有限體積IMEX-RK隱式格式推廣于格子Boltzmann方法，并應用在近空間連續(xù)流區(qū)非規(guī)則航天器解體物繞流模擬，結(jié)論如下：

1）采用適用于近空間連續(xù)流區(qū)殘骸碎片非規(guī)則解體物繞流描述的耦合雙分布函數(shù)格子Boltzmann模型，使其具有普朗特數(shù)、比熱比任意可調(diào)。采用圓函數(shù)替代Maxwell平衡態(tài)分布函數(shù)，發(fā)展了離散速度模型，拓展了格子Boltzmann方法模擬高可壓縮流動問題能力。

2）采用有限體積法數(shù)值離散處理該模型方程以保證構(gòu)造的數(shù)值格式具有更好的通量守恒性及更強的非規(guī)則解體物外形適應能力。為了解決源項剛性問題，選用IMEX-RK格式進行時間項離散。在IMEX-RK格式應用中，利用Boltzmann方程具有的碰撞項對碰撞守恒量的矩為零的特點，推導得到了各階宏觀流動量的顯式表達式，避免了隱式求解各階速度分布函數(shù)。

圖8 CDDF與N-S方法的計算結(jié)果比較Fig.8 Solution com parison between CDDF and N-S solver

3）成功模擬了Riemann問題、雙馬赫反射及翼型繞流問題等案例，驗證了理論模型在可壓流動數(shù)值模擬能力，證實IMEX-RK能很好地解決源項剛性導致的穩(wěn)定性問題。同時看出，所發(fā)展方法在模擬高馬赫數(shù)復雜高溫度梯度流動方面還存在不足，不同的邊界條件對模擬結(jié)果的合理性、數(shù)值穩(wěn)定性會產(chǎn)生影響。本文工作屬初步研究，更精細計算分析驗證有待開展。