高超聲速連續(xù)/稀薄流自適應界面推進重疊網格方法及應用研究

許 嘯，王學德，黃 飛，譚俊杰

(1.南京理工大學能源與動力工程學院，江蘇南京 210094;2.中國航天空氣動力技術研究院，北京 100074)

高超聲速連續(xù)/稀薄流自適應界面推進重疊網格方法及應用研究

許 嘯1，*，王學德1，黃 飛2，譚俊杰1

(1.南京理工大學能源與動力工程學院，江蘇南京 210094;2.中國航天空氣動力技術研究院，北京 100074)

根據高超聲速飛行器在過渡流域飛行時繞流流場具有連續(xù)流與稀薄流并存的特點，提出了一種適用于混合流場自適應界面推進的多尺度重疊網格方法。用兩套獨立生成不同尺度的結構網格和非結構網格覆蓋整個流場。根據連續(xù)/稀薄流界面失效函數(shù)的截止值在結構網格內進行自適應“挖洞”。在挖洞后的稀薄流區(qū)利用非結構網格DSMC方法進行計算;在連續(xù)流區(qū)的結構網格內采用基于N-S方程的CFD算法。兩種算法進行耦合計算和求解，在連續(xù)/稀薄流分界面進行流場信息的插值、交換和處理。根據流場實時耦合計算結果，利用界面失效函數(shù)自適應改變兩套網格的重疊部分。最后，通過對超聲速圓柱和高超聲速帶擴張角圓管繞流的數(shù)值模擬，驗證了方法的可行性和有效性。

高超聲速;連續(xù)/稀薄氣流;自適應界面推進;重疊網格

0 引言

隨著航天技術的不斷發(fā)展，各類飛行器在過渡流域范圍內工作的時間越來越長，其飛行過程中的運動特性逐漸成為研究的熱點。根據Knudsen數(shù)(克努森數(shù))的不同，將飛行器所經歷的流場分為連續(xù)流區(qū)、過渡流區(qū)和自由分子流區(qū)，這其中過渡流區(qū)的求解問題最為復雜。這是因為航天飛行器在過渡流區(qū)飛行時，不僅流場大部分區(qū)域或局部區(qū)域內存在著明顯的稀薄氣體效應，而且由于壓縮效應甚者可能出現(xiàn)連續(xù)流與稀薄流共存的混合流場，從而使整個流場具有連續(xù)與稀薄流混合的特性。過渡流區(qū)飛行器繞流特性問題一直是高超聲氣動技術的難點和熱點。

在過渡流區(qū)，連續(xù)介質模型開始失效，基于N-S方程的CFD(Computational Fluid Dynamics)方法不再適用，而稀薄氣體動力學方法則成為該流域的主要研究方法。這其中Bird提出的DSMC(Direct Simulation of Monte Carlo)方法得到了廣泛的運用［1］。針對連續(xù)與稀薄流動共存的混合流場，則需采用基于連續(xù)流和稀薄流的CFD/DSMC耦合方法進行計算［2-6］。該方法針對不同的區(qū)域采用相應的算法進行計算，既可以完整地反應混合流場的特征，又盡可能地提高了計算效率。在此方面，Hash和Hassan提出了一種CFD/ DSMC的解耦方法，并嘗試采用基于通量的弱耦合方法模擬了Couette流動［7］;T.E.Schwartzentruber和I.D.Boyd對CFD/DSMC耦合方法進行了完善和擴充，以局部Knudsen數(shù)對混合流場進行劃分，采用基于狀態(tài)參數(shù)的方法對高超聲速非平衡氣體流場進行模擬，并提出一種“亞松弛”方法降低統(tǒng)計耗散［8］。李中華等在此基礎上發(fā)展了適用于流場分區(qū)信息交換的亞松弛技術，對CFD/DSMC耦合算法在低克努森數(shù)流場中的應用進行了研究［9］。盡管只采用了單一的網格結構，但這些研究仍為耦合方法的發(fā)展提供了完整的思路。另外在網格方法上，梁杰等以笛卡爾網格DSMC方法為基礎，引入了自適應的非結構網格方法，為重疊網格方法在稀薄氣體中的應用提供了參考［10］。

在混合方法中，由于不同流區(qū)計算方法不同，其分界面的劃定則成為了該方法要解決的首要問題。在連續(xù)/稀薄流混合流場中，很多情況下并沒有明確的界面將流場區(qū)分開來，因此在初始流場中稀薄與連續(xù)區(qū)域的交界面位置都是未知的，需要采用連續(xù)介質失效參數(shù)進行判定。而在流場推進中，取某一時刻的計算結果進行劃定的流場區(qū)域并不一定符合流場最終的情況，因此需要采用自適應方法進行推進，并采用CFD/DSMC耦合方法進行流場計算。連續(xù)流區(qū)采用的CFD方法和稀薄流區(qū)采用的DSMC方法對計算網格的類型、質量、尺度的要求不同，因此采用單一類型的計算網格難以兼顧不同求解方法的特點和要求，也難以得到較為精確可靠的數(shù)值解。此外兩種流動界面的復雜性和推進過程界面的運動性，導致單一的網格類型也難以適應計算。為此，本文提出了一種基于結構和非結構網格的重疊網格界面自適應推進方法，開展了高超聲速連續(xù)/稀薄流混合流場數(shù)值模擬技術研究，編制了通用的計算程序，通過對過渡流域超聲速圓柱繞流和高超聲速帶擴張角圓管繞流的數(shù)值試驗，驗證了算法的可行性與有效性。

1 連續(xù)/稀薄混合流場數(shù)值模擬方法

在連續(xù)流區(qū)域，本文采用基于結構網格N-S方程的CFD算法，以獲得高精度的計算結果和流場分辨率;在稀薄流區(qū)域，采用非結構網格DSMC方法進行計算。

1.1 連續(xù)流區(qū)數(shù)值求解方法

直角坐標系中積分形式的N-S方程為

其中，W為守恒變量，F(xiàn)c為無粘通量，F(xiàn)v粘性通量。對于每一個控制體單元，采用有限體積法對方程進行離散，對流通量離散采用Van Leer矢通量分量格式，并使用Van Albada限制器通過MUSCL方法進行插值。粘性項采用中心差分離散，時間推進采用五步Runge-Kutta法。

連續(xù)流區(qū)CFD方法一般利用正交性較強的結構網格以取得較高精度的解，所用的網格尺度依賴于解算格式和流場的宏觀特性，是基于宏觀層面流場特性的宏觀網格尺度。本文在連續(xù)流區(qū)中采用結構化網格。由于不同流區(qū)之間需要通過單元信息量的插值進行數(shù)值傳遞，因此結構網格單元需要在模擬流場宏觀變化的同時滿足微觀分子運動規(guī)律和統(tǒng)計規(guī)律，使不同區(qū)域間的流場能夠光滑過渡。該網格既需要能夠在插值處模擬出微觀流動變化的情況(尺寸在分子平均自由程左右)，又需要包含有足夠多的分子以滿足稀薄區(qū)域中相關單元的統(tǒng)計需求(至少20個分子)，同時需要能夠模擬近壁面處的粘性效應。據此，本文以來流氣體分子平均自由程的0.1倍至0.33倍為標準設定結構網格單元在近壁面處的尺寸，使其滿足不同區(qū)域間網格的匹配需求。

1.2 稀薄流區(qū)非結構網格DSMC方法

在稀薄流區(qū)采用基于微觀分子氣體動力學的DSMC方法。DSMC方法用有限個模擬分子代替真實氣體分子，在計算機中存儲模擬分子的位置坐標、速度分量以及內能，其值隨模擬分子的運動與邊界的作用以及模擬分子之間的碰撞而不斷變化，通過對計算網格內模擬分子的運動狀態(tài)取樣統(tǒng)計達到求解真實氣體流動的目的。DSMC方法對網格類型依賴程度較小，但對網格尺度要求較高，該方法中網格的作用主要是對流場進行分子微觀取樣和統(tǒng)計。網格尺度過大和過小均得不到正確的計算結果，過大的尺度會導致計算精度的降低，而過小的尺寸會導致統(tǒng)計漲落過于明顯。因此，DSMC方法中計算網格尺度取決于流場的稀薄程度，是基于流場微觀層面上稀薄程度的微尺度網格。一般來說，在生成初始網格時，以來流氣體平均分子自由程的三分之一為標準確定網格單元的大小，這樣可以保證貼近物面的高密度區(qū)的計算精度，同時由于近壁面處分子數(shù)密度大，可以保證該尺寸單元在取樣時有足夠的分子樣本。本文DSMC計算網格采用對復雜外形適應性較強和求解過程易于加密的非結構網格［11-12］，以保證對復雜流域的高度適應性和DSMC方法對網格尺度的要求。

為了解決傳統(tǒng)DSMC方法統(tǒng)計漲落對CFD/ DSMC耦合計算的影響，本文采用基于信息保存IP方法［13］的 DSMC方法。分子模型采用 GSS-3模型［14］，該模型同時具有吸引和排斥的分子作用力，能夠給出正確的氣體輸運性質。利用分子與相鄰單元面的相交理論設計追蹤模擬分子在網格間遷移運動的搜索算法，同時引入碰撞距離的思想來保證正確的模擬結果。

2 自適應推進界面中的連續(xù)失效參數(shù)

連續(xù)與稀薄流混合流場的計算中，分界面的劃分通常采用一個連續(xù)介質失效參數(shù)作為判定依據，并根據實際情況給出截止值，將不同流場部分中連續(xù)失效參數(shù)值與截止值進行比較，大于截止值的部分采用DSMC方法，反之則采用CFD方法，這兩部分的分界面則是CFD/DSMC的交界面。本文采用Boyd提出的局部克努森數(shù)法［15-16］作為界面劃分的依據:

其中，λ是分子平均自由程，Q表示某流動參數(shù)，例如密度，速度或者溫度。通常采用這些局部克努森數(shù)的最大值作為具體計算中的失效參數(shù)，即:

在計算中連續(xù)介質失效參數(shù)的截止值取為0.05。

3 自適應界面推進中的重疊網格方法

3.1 主要思想

本文發(fā)展的連續(xù)/稀薄流混合流場重疊網格法的基本思想是根據不同流區(qū)對計算網格的需要，用兩套獨立生成的不同尺度網格覆蓋整個流場，在不同的流動區(qū)域形成重疊網格。不同流區(qū)的多尺度計算網格一次生成，兩套網格可根據需要實現(xiàn)任意不同類型和尺度網格的重疊、雜交和重疊區(qū)域任意形式的動態(tài)改變。根據連續(xù)/稀薄流界面失效函數(shù)的截止值在連續(xù)流區(qū)進行自適應“挖洞”，動態(tài)改變兩套網格的重疊部分，挖洞后的區(qū)域為DSMC區(qū)域，并激活DSMC區(qū)域的非結構網格利用DSMC方法進行計算。同時根據自適應過程，對于不再參與稀薄流區(qū)計算的非結構網格單元進行“休眠”，并激活相應位置的結構網格單元參與連續(xù)流區(qū)的計算。在連續(xù)/稀薄流區(qū)分界面附近的重疊區(qū)域利用插值方法實現(xiàn)流場信息的傳遞。

3.2 重疊網格技術中的區(qū)域挖洞方法

由于結構網格和非結構網格單元均是獨立覆蓋全流場的，二者數(shù)據結構完全不同，當根據失效參數(shù)的截止值以結構網格進行區(qū)域洞邊界劃分時，如何確定非結構網格的點是否落在洞內成為了該方法的關鍵。

由于非結構網格的拓撲結構的無序性，本文以結構網格為背景網格發(fā)展了一種區(qū)域挖洞方法，其一般過程為:

(1)首先以結構網格作為背景網格，遍歷所有非結構網格單元，利用其形心坐標找出并記錄所屬的結構網格單元，將結構網格單元和非結構網格單元建立起對應關系，即確定結構網格單元與非結構網格單元的包含關系。

(2)當進行重疊網格洞邊界劃分時，以流場計算的結構網格為基礎進行洞邊界劃分，通過網格點上的局部克努森數(shù)的大小將點分為計算點(激活狀態(tài))與非計算點(休眠狀態(tài))，將與非計算點相鄰的計算點設置為邊界點，根據單元形心點的屬性將單元分為計算單元與非計算單元，并從計算單元的邊界處起，感染相鄰的非計算單元作為插值單元。此時初始洞邊界已完成。

(3)遍歷結構網格所有洞內(DSMC區(qū)域)和洞外單元(CFD區(qū)域)，分別將所有包含在內的非結構網格點相應的設置為計算點(激活狀態(tài))和非計算點(休眠)。再檢查非結構網格中的所有單元，根據單元中所有點的具體屬性(計算點、非計算點和界面點)將其設置為洞內單元，洞外單元和插值單元。

(4)遍歷所有網格單元，根據結構和非結構網格單元屬性給出CFD和DSMC計算區(qū)域，如圖1所示。

圖1 挖洞后CFD/DSMC計算網格單元Fig.1 Overlapping grid of CFD/DSMC after hole cutting

3.3 重疊網格插值單元的信息交換方法

在計算中，DSMC和CFD區(qū)域的信息傳遞通過插值單元處的信息轉換完成。由于DSMC的計算結果為統(tǒng)計量，在信息轉換時不可避免的存在統(tǒng)計耗散，本文采用IP方法對其進行改進，有效地抑制了DSMC統(tǒng)計漲落對CFD的影響。根據CFD/DSMC方法交界面兩側的信息傳遞方式，可將耦合方法分為兩種:基于通量法和基于狀態(tài)參數(shù)法?；谕糠ㄒ蕾囉诖罅康姆肿咏y(tǒng)計，在分子數(shù)不足的情況下可能會發(fā)生計算崩潰的情況，因此本文采用基于狀態(tài)參數(shù)法進行計算。

狀態(tài)參數(shù)法首先對洞邊界兩側不同區(qū)域單元內的狀態(tài)參數(shù)進行轉換，生成洞邊界外側單元的信息值，然后使用該值進行流場計算。但在結構與非結構網格互相重疊的情況下，不同區(qū)域的洞邊界面并不重合，因此首先需要根據各區(qū)域的洞邊界位置確定其外側單元，然后以該單元在相鄰區(qū)域中的貢獻單元信息為基礎進行插值，最后將計算出的單元信息轉換為洞邊界上各計算方法所需要的物理量。這樣使得各區(qū)域在處理洞邊界的通量信息時，外側單元與計算單元之間的界面相互重合，從而保證了該區(qū)域內的通量守恒。

當CFD區(qū)域為DSMC區(qū)域提供邊界條件時，即重疊區(qū)數(shù)值需要進行從宏觀量到微觀量的轉換，此時連續(xù)流區(qū)的宏觀參數(shù)為信息源，其信息交換過程如下:

(1)確定在連續(xù)流重疊區(qū)的非結構網格單元，標記需要進行信息重置的網格。

(2)對所有被標記到的非結構網格單元進行信息重置，將該單元內分子的速度，溫度，位置等信息歸零，刪除分子鏈表中相關信息。

(3)對重疊區(qū)中相關結構網格單元進行信息轉換處理，即以結構網格單元內的密度，速度，溫度等宏觀信息量為基礎，用Chapman-Enskog理論方法生成結構網格單元內的分子信息。

(4)在結構網格單元內對分子位置進行隨機分布，將生成的新分子納入到所對應非結構網格單元的分子信息鏈表中。

重疊區(qū)內當DSMC區(qū)域向CFD區(qū)域提供邊界條件時，即重疊區(qū)數(shù)值需要進行從微觀量到宏觀量的轉換，此時稀薄流區(qū)的微觀參數(shù)為信息源，信息交換過程如下:

(1)對重疊區(qū)中相關結構網格單元內的宏觀量進行重置，并根據單元間的關系鏈表確定與其相關的非結構網格單元編號。

(2)對每個結構網格單元，首先統(tǒng)計該單元的所有相關非結構網格單元。然后分別求出各相關非結構網格單元中心到該單元中心的距離，將該值的倒數(shù)作為加權系數(shù)，按照如下公式計算:

式中，L為距離的倒數(shù)，Qi為信息量。求出Q值之后，將其作為插值單元中的信息值進行通量計算。

(3)按上述方法整理所有結構網格單元中新生成的宏觀量，將其賦值給相關的結構網格單元。

4 連續(xù)/稀薄流自適應界面算法流程

本文所設計的連續(xù)/稀薄流自適應界面推進重疊網格算法的主要流程如下:

(1)采用結構網格和非結構網格分別布置全流場，并采用結構網格CFD方法進行流場初始計算，以迭代Nstep步后的流場作為耦合開始計算的初始流場。

(2)采用局部克努森數(shù)作為稀薄與連續(xù)流區(qū)域劃分的標準，其截止值取0.05，在結構網格內劃出初始的洞邊界。

(3)以結構網格初始洞邊界為基準，在非結構網格內進行挖洞，劃出相應的DSMC計算、插值和非計算區(qū)域。

(4)進行CFD和DSMC的耦合迭代計算，由于耦合迭代步數(shù)較小時界面移動不明顯，為此本文在程序中設定自適應的循環(huán)迭代間隔數(shù)NADAP。

(5)如果MOD(NPR，NADAP)=0進入自適應階段，按照新計算的流場結果再次劃定結構和非結構網格的洞邊界區(qū)域，這里NPR為循環(huán)總數(shù)，NADAP自適應循環(huán)間隔數(shù)。

(6)統(tǒng)計結構網格與非結構網格的洞邊界情況，當上下兩個自適應間隔時間步內的洞邊界位置不再變化時，關閉自適應挖洞的過程，將該分界面作為最終分界面。按此界面耦合計算至收斂，否則返回(4)繼續(xù)進行流場的自適應界面推進計算。

(7)由于重疊網格的不同區(qū)域算法不同，因此其判定標準也不相同。NS方法一般以密度作為殘差的標準，計算時首先記錄所有單元兩次時間步之間的密度相差值之和，將每一步的殘差和與初始殘差和相比，當該比值小于或者等于10-4時認為進入收斂狀態(tài);DSMC方法一般以統(tǒng)計漲落的情況作為判定收斂的標準，首先統(tǒng)計兩個統(tǒng)計循環(huán)中分子總數(shù)的相差值，當該值與分子總數(shù)的統(tǒng)計平均值之比小于10-4時，認為計算已經收斂。在重疊網格耦合計算中，當這兩個條件均滿足時，認為流場計算收斂。

5 數(shù)值算例與分析

5.1 超聲速圓柱繞流

本文首先計算了超聲速圓柱繞流來驗證算法的有效性，計算模型來自文獻［17］，圓柱直徑為0.1 m，采用結構和非結構網格覆蓋全流場，非結構網格共由41 694個三角形單元和21090個點組成;結構網格規(guī)模為50×100個計算單元。計算條件同文獻［17］，來流氣體為純氮氣，粘性系數(shù)的溫度冪指數(shù)為0.74，分子自由度為5，來流馬赫數(shù)為4，速度為1 412.5 m/s，溫度為300 K，圓柱壁面溫度為500 K，密度為4.65× 10-6kg/m3。流場中分子數(shù)密度為1.0×1020，每個模擬分子所代表的真實分子數(shù)設為2.5×1013，時間步長為4.0×10-6s。

圖2給出了界面自適應過程中連續(xù)/稀薄流區(qū)計算網格單元變化情況，從圖中可以看出隨著流場的發(fā)展和激波的形成及移動，連續(xù)流區(qū)與稀薄流區(qū)逐漸發(fā)展的過程。在弓形激波和近壁面附近由于流場物理量梯度的急劇變化導致稀薄氣體效應明顯，生成了局部范圍內的DSMC區(qū)域，并且隨著弓形激波的發(fā)展，激波附近的DSMC區(qū)域逐漸擴展。在圓柱的尾流區(qū)，由于低密度效應生成了局部DSMC區(qū)域。從圖中可以看出，尾流區(qū)的DSMC區(qū)域呈現(xiàn)先增大后減小的趨勢，這主要是因為在流場初始計算時，尾部區(qū)域流場速度取為來流均勻速度，較大的初始速度造成計算過程中尾部區(qū)域首先出現(xiàn)大面積的稀薄流區(qū)，而后由于氣流對尾部稀薄流域的氣體補充作用，稀薄流區(qū)又逐漸減小。

圖2 圓柱繞流連續(xù)/稀薄流區(qū)計算網格自適應界面推進過程Fig.2 Overlapping grids during adaptive interface advancement process of a cylinder supersonic flow

圖3和圖4分別給出了圓柱繞流自適應界面推進過程中密度和溫度的計算等值線圖，并分別指出了等值線和不同區(qū)域分界面的位置。從圖中可以看出，在流場的發(fā)展過程中，整個流場的分辨率較高，圓柱頭部弓形激波發(fā)展過程中激波面清晰，這表明本文采用基于信息保存方法的DSMC方法有效抑制了統(tǒng)計漲落對CFD/DSMC耦合過程的影響。同時也表明在自適應推進過程中本文所使用的連續(xù)流CFD方法與稀薄流DSMC方法耦合邊界處信息交換方法是可行的。

圖3 圓柱繞流連續(xù)/稀薄流區(qū)自適應界面推進過程密度計算等值線圖Fig.3 Density contours of supersonic cylinder during adaptive interface advancement process

圖4 圓柱繞流連續(xù)/稀薄流區(qū)自適應界面推進過程溫度計算等值線圖Fig.4 Temperature contours of supersonic cylinder during adaptive interface advancement process

圖5給出了圓柱繞流壁面氣動系數(shù)計算結果比較，圖中橫坐標L為以駐點為起始點的弧長值，從圖中可以看出本文計算結果與文獻［17］以及DSMC計算結果吻合得較好，進一步表明本文發(fā)展的方法是可行的。圖5(a)為沿上壁面的壓力系數(shù)分布，可以看出從駐點起壓力系數(shù)分布類似于余弦函數(shù)曲線分布，與理論分析結果一致，圖5(b)給出從駐點沿上壁面的摩擦系數(shù)分布的計算比較，在駐點處摩擦系數(shù)為零，沿上壁面逐漸增加，在圓心角45°左右增加至最高峰后沿圓柱壁面又減小，與理論分析結果一致。

圖5 圓柱壁面物理量分布計算結果比較Fig.5 Comparison of calculated results on the cylinder wall

5.2 帶有擴張角的高超聲速圓管繞流

帶擴張角的圓管結構示意如圖6所示，計算區(qū)域為管外的流場，即圓管外壁與擴張角形成的組合體。采用結構和非結構網格覆蓋全流場，非結構網格共由76 407個三角形單元和39406個點組成;結構網格為由60×120個計算單元組成。

圖6 帶有擴張角的圓管結構示意圖Fig.6 The structure of the hollow cylinder-flare

計算條件同文獻［18］，來流氣體為純氮氣，粘性系數(shù)的溫度冪指數(shù)為0.74，分子自由度為5，來流馬赫數(shù)為10.5，速度為2 301.7 m/s，溫度為118.2 K，密度為9.023×10-4kg/m3。分子數(shù)密度為1.94×1022，每個模擬分子所代表的真實分子數(shù)設為3.5×1014，時間步長為2.0×10-6s。

圖7 自適應過程完成后計算區(qū)域網格劃分及密度云圖Fig.7 Overlapping grid and density contours of after adaptive interface advancing process

圖7給出了自適應界面推進完成后連續(xù)流區(qū)與稀薄流區(qū)計算網格及密度云圖。從圖7(b)中可以看出，圓管壁面斜激波和擴張壁面斜激波的相互干擾特征明顯，激波與激波同側相交處密度增高，整個繞流顯示了較高的流場分辨率，進一步說明了本文所設計的自適應方法是可行的。

圖8給出了擴張圓管摩擦系數(shù)與壁面熱流分布計算結果的比較，從圖中可以看出本文摩擦系數(shù)和熱流計算結果與文獻［18］中給出的計算值與實驗值吻合得較好。

為了考察自適應界面推進重疊網格法的計算效率，本文將算例1與算例2在相同配置的計算機上分別利用非結構網格DSMC方法和本文方法進行了計算效率比較。所用計算機為曙光CB85-F四路四核高性能刀片計算服務器，CPU為 AMD Opteron 8374HE，主頻2.2 GHz，16 GB內存。表1列出了每種方法從計算開始至收斂所用計算時間的比較，表2列出了統(tǒng)計穩(wěn)定后至收斂時的計算時間和各區(qū)域在總流場區(qū)域中所占的比例。從中可以看出與DSMC方法相比較，本文方法有效地減少了稀薄區(qū)域的計算量，降低了計算時間，提高了計算效率。

圖8 擴張圓管壁面物理量分布計算結果比較Fig.8 Comparison of calculated results on the wall of hollow cylinder-flare

表1 總計算時間的比較Table 1 Comparison of whole computing time

表2 統(tǒng)計后至收斂的計算時間和各區(qū)域所占比例Table 2 The calculation time from statistics to convergence and the proportion of different regions of the grid

為了比較兩種算法的收斂歷程，圖9給出了算例1中圓柱的阻力系數(shù)收斂變化曲線。從圖中可以看出，與DSMC方法相比，本文方法不僅在計算速度上有明顯的優(yōu)勢，并且阻力系數(shù)變化和收斂過程更加穩(wěn)定，進一步表明了混合方法的優(yōu)勢。

圖9 圓柱表面阻力系數(shù)的收斂過程Fig.9 The convergent process of the drag coefficients on the cylindrical surface

6 結論

本文開展了過渡流區(qū)高超聲速連續(xù)/稀薄流混合流場數(shù)值模擬技術研究。提出了一種基于結構和非結構網格基礎上的重疊網格界面自適應推進方法。通過超聲速圓柱繞流和高超聲速擴張圓管繞流的數(shù)值模擬驗證，可以得到如下結論:

(1)數(shù)值計算結果表明該方法在處理連續(xù)/稀薄流混合流場的模擬問題是可行的，所得流場分辨率高，計算結果與文獻和實驗比較吻合較好。

(2)該方法在求解過程中可以自適應捕捉連續(xù)流與稀薄流的分界面，實現(xiàn)了界面捕捉和耦合計算過程的統(tǒng)一。

(3)通過計算效率的比較，該方法與非結構網格DSMC方法相比，能夠有效降低計算時間，提高計算效率。

［1］ Bird G A.Molecular gas dynamics and the direct simulation of gas flows［M］.Oxford:Clarendon Press，1994.

［2］ Thomas E S，Leonardo C S，et al.Modular implementation of a hybrid DSMC-NS algorithm for hypersonic non-equilibrium flows［R］.AIAA 2007-613.

［3］ George J D，Boyd I D.Simulation of nozzle plume flows using a combined CFD/DSMC approach［R］.AIAA-99-3454，1999.

［4］ Wadsworth D C，Erwin D A.Two-dimensional hybrid contiuum particle simulation approach for rarefied flows［R］.AIAA-92-2975，1992.

［5］ Schwartzentruber T E，Scalabrin L C，Boyd I D.Hybrid particle continuum simulations of non-equilibrium hypersonic blunt body flow fields［R］.AIAA 2006-3602.

［6］ Li Zhihui，Li Zhonghua，et al.Coupled simulation of mixture plume for attitude control satellite thruster［J］.Acta Aerodynamica Sinica，2012，30(4):483-491.(in Chinese)

李志輝，李中華，等.衛(wèi)星自控發(fā)動機混合物羽流場分區(qū)耦合計算研究［J］.空氣動力學學報，2012，30(4):483-491.

［7］ Hash D B，Hassan H A.Assessment of schemes for coupling Monte Carlo and Navier-Stokes solution methods［J］.J.Thermophys Heat Transf.，1996，10:242-249.

［8］ Schwartzentruber T E，Boyd I D.A hybrid particle-continnum method applied to shock waves［J］.Journal of Computational Physics，2006，215:402-416.

［9］ Li Zhonghua，Li Zhihui，Li Haiyan，et al.Research on CFD/ DSMC hybrid numerical method in rarefied flows［J］.Acta Aerodynamica Sinica，2013，2013，31(3):282-287.(in Chinese)

李中華，李志輝，李海燕，等.過渡流區(qū)N-S/DSMC耦合計算研究［J］.空氣動力學學報，2013，31(3):282-287.

［10］Liang Jie，Yan Chao，Du Boqiang.An algorithm study of three-dimensional DSMC simulation based on two-level cartesian coordinates grid structure［J］.Acta Aerodynamica Sinica，2010，28(4):466-471.(in Chinese)

梁杰，閻超，杜波強.基于兩級直角網格結構的三維DSMC算法研究［J］.空氣動力學學報，2010，28(4):466-471.

［11］Wang Xuede，Wu Yizhao，Xia jian.Strategies for implement-tation of 2D unstructured DSMC method and its application［J］.Acta Aerodynamica Sinica，2007，25(1):110-115.(in Chinese)

王學德，伍貽兆，夏健.一類二維非結構網格DSMC方法實現(xiàn)策略及其應用［J］.空氣動力學學報，2007，25(1):110-115.

［12］Wang Xuede，Wu Yizhao，Xia jian.A parallel algorithm of 2D Unstructured DSMC method with dynamic load balance［J］.Acta Aerodynamica Sinica，2007，25(3):339-344.(in Chinese)

王學德，伍貽兆，夏健.動態(tài)負載平衡的二維非結構網格DSMC并行算法研究［J］.空氣動力學學報，2007，25(3):339-344.

［13］Sun Q H，Boyd I D.Theoretical development of the information preservation method for strongly nonequilibrium gas flows［C］//The 38th AIAA Thermophysics Conference，Toronto，Ontario，Canada，June 2005.

［14］ Wang Xuede，Wu Yizhao，Xia jian.GSS-3 molecular model in Monte Carlo simulation［J］.Chinese Journal of Computational Physics，2009，26(1):64-70.(in Chinese)

王學德，伍貽兆，夏健.蒙特卡羅方法中分子作用模型的設計及應用［J］.計算物理，2009，26(1):64-70.

［15］Wang W L，Boyd I D.Continuum Breakdown in Hypersonic Viscous Flows［R］.AIAA 2002-0651.

［16］Lofthouse A J，Boyd I D，Wright M J.Effects of continuum breakdown on hypersonic aerothermodynamics［R］.AIAA-06-0993，2006.

［17］王學德.高超聲速稀薄氣流非結構網格DSMC及并行算法研究［D］.南京:南京航空航天大學，2006.

［18］Wang W L，Boyd I D.Hybrid DSMC-CFD simulations of hypersonic flow over sharp and blunted bodies［R］.AIAA 2003-3644.

An overlapping grid adaptive interface advancement method for hypersonic continuum-rarefied flow and its application

Xu Xiao1，*，Wang Xuede1，Huang Fei2，Tang Junjie1
(1.School of Energy and Power Engineering，NUST，Nanjing 210094，China; 2.China Academy of Aerospace Aerodynamics Technology，Beijing 100074，China)

When the hypersonic aircraft flies in the transitional regime，both continuum and rarefied flows appear in the flow field around the aircraft.A multi-scale adaptive advanced overset grid method is proposed to deal with such a flow problem.The whole flow field is covered by the structured and unstructured grids which are independently generated in different scales.The adaptive hole cutting process is advanced by the cut-off value of the interface breakdown function of the continuum and rarefied flows in the structured grid.When the hole cutting process is completed，the rarefied regions are computed using the unstructured DSMC method while the continuum regions are computed using the structured CFD method based on the N-S equations.To couple the two algorithms，the flow information at the boundary between the continuum and rarefied regions is interpolated，exchanged and processed.According to the instantaneous solution of the coupled flow field，the overlapping parts of the different grids are advanced by the interface breakdown function.The feasibility and effectiveness of the present method are validated by the simulation results of the supersonic flow across a cylinder and hypersonic flow across a hollow cylinderflare.

hypersonic flow;continuum-rarefied flow;adaptive interface advancement;overlapping grid

V211.3

Adoi:10.7638/kqdlxxb-2014.0015

0258-1825(2015)04-0563-09

2014-03-14;

2014-07-18

南京理工大學紫金之星資助計劃;江蘇省青藍工程資助項目

許嘯*(1986-)，男，江蘇興化人，博士研究生，主要從事稀薄氣體動力學研究.E-mail:xuxiao805@163.com

許嘯，王學德，黃飛，等.高超聲速連續(xù)/稀薄流自適應界面推進重疊網格方法及應用研究［J］.空氣動力學學報，2015，33(4): 563-571.

10.7638/kqdlxxb-2014.0015 Xu X，Wang X D，Huang F，et al.An overlapping grid adaptive interface advancement method for hypersonic continuum-rarefied flow and its application［J］.Acta Aerodynamica Sinica，2015，33(4):563-571.