那旭東，夏智勛，馬立坤，顏小婷

(國防科技大學(xué) 空天科學(xué)學(xué)院，長沙 410073)

0 引言

固體火箭發(fā)動機(jī)是火箭的動力裝置，由于其具有結(jié)構(gòu)簡單、使用方便、可靠性高、成本低等特點(diǎn)，常作為一級或多級的動力裝置被廣泛應(yīng)用于運(yùn)載火箭、戰(zhàn)術(shù)武器以及戰(zhàn)略武器等方面[1]。固體火箭發(fā)動機(jī)的設(shè)計(jì)涉及推進(jìn)劑的點(diǎn)火與燃燒，推進(jìn)劑、殼體、絕熱層及噴管的結(jié)構(gòu)力學(xué)，內(nèi)部流場及外部羽流的流體動力學(xué)，含能材料物理化學(xué)特性，組件的老化損壞以及不同潛在失效模型分析等一系列學(xué)科[2]。此外，固體火箭發(fā)動機(jī)工作過程涉及不同空間(μm～m)、時(shí)間(10-6～102s)尺度的燃燒及流動[3]，在燃燒室內(nèi)點(diǎn)火器點(diǎn)火產(chǎn)生出高溫高壓的氣體點(diǎn)燃固體推進(jìn)劑，推進(jìn)劑燃燒的氣體產(chǎn)物生成壓力作用于固體藥柱表面，促使藥柱變形，另外氣體產(chǎn)物向藥柱傳遞熱量，促進(jìn)推進(jìn)劑的快速燃燒。同時(shí)，藥柱的變形與快速燃燒又影響著燃燒室內(nèi)氣體流動與壓力分布等[4]。因此，固體火箭發(fā)動機(jī)的工作過程是一個(gè)復(fù)雜的、強(qiáng)耦合的物理化學(xué)過程。

自20世紀(jì)70年代以來，采用計(jì)算機(jī)仿真技術(shù)對固體火箭發(fā)動機(jī)進(jìn)行輔助設(shè)計(jì)及工作過程模擬已成為固體火箭發(fā)動機(jī)研究的一個(gè)重要手段。近年來，隨著計(jì)算機(jī)性能的提升以及CFD技術(shù)的發(fā)展，各國對固體火箭發(fā)動機(jī)工作過程的數(shù)值模擬也變得越來越精細(xì)，從零維到多維、從不可壓到可壓、從定常到非定常、從無粘到有粘、從單相到多相等各個(gè)方面不斷發(fā)展[5]。各國在對固體火箭發(fā)動機(jī)工作過程進(jìn)行數(shù)值仿真時(shí)，除了采用Fluent、Star-CD及ABAQUS等商業(yè)軟件外，還開發(fā)了許多相應(yīng)的設(shè)計(jì)仿真程序和軟件[6-12]，并取得了一系列可喜的成果，其中以美國先進(jìn)火箭仿真中心CSAR(Center for Simulation of Advanced Rockets)開發(fā)的Rocstar軟件包最具有代表性[13]。

Rocstar能夠?qū)腆w火箭發(fā)動機(jī)從點(diǎn)火到穩(wěn)態(tài)燃燒工作過程進(jìn)行三維、全耦合及高精度仿真模擬，求解還包括鋁顆粒、煙塵和各種化學(xué)成分的流場[14-15]。本文主要就CSAR的由來，Rocstar的發(fā)展歷程、程序架構(gòu)組成、求解模塊及對應(yīng)的算法、程序驗(yàn)證和求解流程幾方面對Rocstar進(jìn)行介紹。

1 Rocstar開發(fā)背景、機(jī)構(gòu)組成及發(fā)展歷程

1.1 開發(fā)背景和機(jī)構(gòu)組成

20世紀(jì)90年代，美國為了保持戰(zhàn)略平衡，隸屬于能源部的國家核安全管理局發(fā)起了先進(jìn)模擬與計(jì)算ASC (Advanced Simulation and Computing Program)項(xiàng)目[16-17]。ASC項(xiàng)目的目標(biāo)是提升先進(jìn)的計(jì)算建模和仿真能力來實(shí)現(xiàn)對復(fù)雜、多組件系統(tǒng)的大規(guī)模并行仿真模擬，以便在沒有地下核試驗(yàn)的情況下，通過數(shù)值模擬試驗(yàn)來檢測其核武器的安全性與可靠性。

伊利諾斯大學(xué)為響應(yīng)ASC計(jì)劃，于1997年成立了先進(jìn)火箭仿真中心CSAR[18]。CSAR希望開發(fā)一套固體火箭發(fā)動機(jī)虛擬樣機(jī)工具，用于開展固體火箭發(fā)動機(jī)主要組成部件及部件間動態(tài)相互作用的詳細(xì)建模和仿真，從而對整個(gè)固體火箭發(fā)動機(jī)工作狀態(tài)進(jìn)行精確的預(yù)測[19-20]。

為完成固體火箭發(fā)動機(jī)全系統(tǒng)仿真的這一終極目標(biāo)，CASR將開發(fā)人員分成了燃燒和含能材料組、計(jì)算科學(xué)組、流體力學(xué)組、結(jié)構(gòu)和材料組4個(gè)大的研究組以及為了完成特定問題而建立的其他小組。從1997年創(chuàng)建CSAR到2006年底，將近8年時(shí)間，CSAR完成了立項(xiàng)之時(shí)所設(shè)定的目標(biāo)[4]。此后，Illinois Rocstar LLC作為一個(gè)獨(dú)立的商業(yè)部門和IR&D公司開始負(fù)責(zé)Rocstar的開發(fā)和技術(shù)支持，但仍由美國能源部提供資金支持[13,21]。

1.2 發(fā)展歷程及可實(shí)現(xiàn)的功能

CSAR在開發(fā)Rocstar程序過程中采用了逐步深入的開發(fā)策略[22]，圖1給出了Rocstar各個(gè)版本所能處理的幾何及物理模型的復(fù)雜程度。隨著版本的不斷提高，其所能處理物理問題也越來越精細(xì)，幾何模型也越來越復(fù)雜。

當(dāng)版本發(fā)展到Rocstar 3時(shí)，程序中已包含有許多物理應(yīng)用求解器，這些求解器通過一個(gè)強(qiáng)大的耦合框架整合在一起。所有求解器都可進(jìn)行高效且大規(guī)模的并行求解，從而可以對異常復(fù)雜的問題進(jìn)行精細(xì)化的仿真計(jì)算[23]。圖2為Rocstar 3中一系列可供用戶選擇的計(jì)算模型和算法。

圖1 Rocstar各版本功能Fig.1 Functions of Rocstar in each version

圖2 Rocstar 3中可用的模型和算法Fig.2 Models and algorithms in Rocstar 3

在Illinois Rocstar LLC接手Rocstar 3后，將其最新一代的多物理場仿真軟件包Rocstar 3代碼開源[24-25]，并將其命令為Rocstar仿真套件(Rocstar Simulation Suite)。

2 Rocstar架構(gòu)及模塊功能

2.1 Rocstar開發(fā)思想及架構(gòu)組成

固體火箭發(fā)動機(jī)是一種性能優(yōu)越的火箭動力裝置，其結(jié)構(gòu)簡單，其主要由點(diǎn)火器、殼體、推進(jìn)劑、絕熱層和噴管等部件組成，如圖3所示。

固體火箭發(fā)動機(jī)的工作本質(zhì)是燃燒流動過程進(jìn)行的兩次能量轉(zhuǎn)換，其機(jī)理復(fù)雜、環(huán)節(jié)眾多：如在燃燒室多相流動過程中就涉及金屬分散燃燒、凝相破碎/聚合、熔渣沉積、熱結(jié)構(gòu)燒蝕、多相輻射對流熱傳導(dǎo)等過程，如圖3所示。此外，組成發(fā)動機(jī)各部件間的耦合是強(qiáng)非線性的，工作過程涉及不同空間、時(shí)間尺度的燃燒及流動，因此其工作過程是一個(gè)復(fù)雜的、強(qiáng)耦合的物理化學(xué)過程，涉及大量的學(xué)科和科學(xué)問題。

圖3 固體火箭發(fā)動機(jī)基本組成部件及相關(guān)物理化學(xué)過程Fig.3 Basic components and physicochemical processes of SRM

從基本原理出發(fā)，進(jìn)行詳細(xì)、全部件耦合的固體火箭發(fā)動機(jī)仿真模擬，需要開發(fā)代碼來實(shí)現(xiàn)每個(gè)火箭部件的數(shù)學(xué)模型以及它們之間的相互作用。這不僅需要更大的計(jì)算能力，而且要求程序中各求解模塊的數(shù)值正確，且相互保持一致的物理作用。因此，進(jìn)行固體火箭發(fā)動機(jī)全系統(tǒng)的建模和仿真極為復(fù)雜。

CSAR在開發(fā)Rocstar時(shí)，多物理場系統(tǒng)集成的方法是開發(fā)一個(gè)主程序和一套接口代碼，以及一系列可單獨(dú)運(yùn)行的求解程序，然后通過接口代碼實(shí)現(xiàn)主程序?qū)ο鄳?yīng)求解程序的調(diào)用[4]。Rocstar仿真應(yīng)用程序的體系結(jié)構(gòu)見圖4，包括了大量用于對流體-固體-燃燒三者之間相互作用、耦合的組件程序。

圖4 Rocstar仿真套件體系結(jié)構(gòu)Fig.4 Rocstar simulation suite architecture

從圖4可看出，Rocstar主要分為前處理、物理模塊求解器、集成架構(gòu)以及后處理4個(gè)功能區(qū)域，其中集成框架是整個(gè)Rocstar的核心所在，包含了許多用于支持多物理場耦合仿真的服務(wù)及協(xié)調(diào)模塊。

在Rocstar開發(fā)過程中，采用了面向?qū)ο蟮某绦蜷_發(fā)思想，服務(wù)及協(xié)調(diào)模塊采用C++語言編寫，物理應(yīng)用求解器采用Fortran 90語言編寫并采用MPI實(shí)現(xiàn)并行仿真。Rocstar可在集群平臺上進(jìn)行固體火箭發(fā)動機(jī)一體化、全系統(tǒng)耦合的大規(guī)模仿真計(jì)算[22]。

2.2 協(xié)調(diào)及服務(wù)模塊

Roccom 是Rocstar最底層的框架。Roccom可實(shí)現(xiàn)不同程序模塊中數(shù)據(jù)及函數(shù)交互、內(nèi)部程序的輸入輸出(Rocin,Rocout)、網(wǎng)格之間的通信操作(Rocsurf,Rocmap)等功能，其中程序模塊可支持使用不同程序語言，主要為C++和F90。通過調(diào)用Roccom，物理模塊相互之間可以進(jìn)行通信，同時(shí)也可以訪問一系列服務(wù)，包括網(wǎng)格修改、數(shù)據(jù)傳遞、數(shù)據(jù)輸入輸出、性能分析及負(fù)載平衡等。

Rocman是Rocstar的主要驅(qū)動和協(xié)調(diào)服務(wù)模塊，主要執(zhí)行和管理不同組件程序間的相關(guān)作用，控制物理應(yīng)用的執(zhí)行過程，包含數(shù)值仿真的初始化、耦合時(shí)間步、交界面處的跳躍條件以及輸出存儲和仿真終止等。Rocman最新的版本為Rocman 3，采用C++編寫，其實(shí)現(xiàn)的功能比以前更加完善和通用。

Rocprop、Rocmop、Rocrem及Rocface用于耦合界面演化與數(shù)據(jù)傳遞。

Rocprop是不同域之間接觸面演化服務(wù)模塊，主要為Rocstar提供拉格朗日面追蹤及演化服務(wù)，用于計(jì)算推進(jìn)劑表面由于燃燒所帶來的退移過程。Rocprop既可以用來進(jìn)行多域耦合仿真也可以單獨(dú)進(jìn)行流體仿真或結(jié)構(gòu)力學(xué)仿真。Rocprop 實(shí)現(xiàn)了兩種分界面演化算法：(1)早期采用的粒子標(biāo)記方法；(2)更先進(jìn)的面位移算法。面位移算法(FOM)對于處理邊角附近的面的運(yùn)動時(shí)很有優(yōu)勢[26]。

對于仿真過程涉及諸如固體藥柱燃面退移、藥柱裂紋動態(tài)發(fā)展等幾何構(gòu)型發(fā)生改變的情況，為保證網(wǎng)格的質(zhì)量，在Rocstar中有兩種網(wǎng)格修改方案來處理不同等級的網(wǎng)格失效情況。非結(jié)構(gòu)網(wǎng)格的光順(不改變網(wǎng)格節(jié)點(diǎn)數(shù)量及節(jié)點(diǎn)連接情況)是通過Rocmop調(diào)用Mesquite程序包完成的，Mesquite[27]是由美國桑迪亞國家實(shí)驗(yàn)室開發(fā)的網(wǎng)格算法程序。每一個(gè)分區(qū)同時(shí)調(diào)用Mesquite程序，用于實(shí)時(shí)提供真實(shí)和虛擬(在計(jì)算域外)的節(jié)點(diǎn)。如果光順不能提高網(wǎng)格的質(zhì)量，則采用Rocrem模塊對網(wǎng)格進(jìn)行局部的修改或重新進(jìn)行網(wǎng)格生成，Rocrem通過調(diào)用商業(yè)軟件包Simmetrix來完成這一工作。當(dāng)網(wǎng)格重新生成后，采用Rocface模塊將計(jì)算結(jié)果映射到新的網(wǎng)格上繼續(xù)進(jìn)行計(jì)算。Rocface[28]是一個(gè)數(shù)據(jù)映射服務(wù)模塊，可將獨(dú)立離散面之間的求解數(shù)據(jù)進(jìn)行傳送。

2.3 Rocstar前處理

對于幾何建模，只要CAD軟件可生成被網(wǎng)格劃分軟件讀取的幾何模型就可以采用。CSAR開發(fā)人員采用Pro/Engineer來生成流體及固體域的CAD模型，并輸出為通用幾何模型接口IGES格式。對于網(wǎng)格劃分，CSAR開發(fā)人員指定了Rocstar前處理中的網(wǎng)格及邊界條件的格式。目前，對于流體力學(xué)程序采用Gridgen軟件來生成網(wǎng)格及邊界條件，對于結(jié)構(gòu)力學(xué)程序采用Partran或Truegrid來生成對應(yīng)的網(wǎng)格及邊界條件。當(dāng)網(wǎng)格和邊界條件以指定格式輸出后，在輸入數(shù)據(jù)集準(zhǔn)備工具Rocprep的輔助下為每個(gè)物理仿真應(yīng)用創(chuàng)建數(shù)據(jù)輸入集以及為并行運(yùn)算進(jìn)行分區(qū)劃分。

2.4 求解器

2.4.1 流體力學(xué)求解器

Rocstar中流體動力學(xué)求解器主要有兩個(gè)，分別為Rocflo和Rocflu，主要用于求解可壓縮多相流。兩個(gè)求解器均采用基于體心的有限體積離散方法求解。

Rocflo采用任意拉格朗日-歐拉ALE (Arbitrary Lagrangian Eulerian)算法在多塊結(jié)構(gòu)網(wǎng)格上進(jìn)行求解[28]，ALE可以考慮由于推進(jìn)劑變形或燃燒帶來的邊界運(yùn)動情況。Rocflo求解器空間離散格式采用中心格式或二階TVD格式。時(shí)間離散，除了顯示的Runge-Kutta方法，Rocflo還可以使用一個(gè)雙時(shí)間步算法來處理比CFL條件限制下的時(shí)間步長更大的時(shí)間步長。Rocflo的并行劃分是基于塊拓?fù)浣Y(jié)構(gòu)，每個(gè)處理器包含1個(gè)或多個(gè)網(wǎng)格塊。

Rocflu可在非結(jié)構(gòu)四面體網(wǎng)格或混合四面體/六面體/金字塔/棱柱網(wǎng)格上進(jìn)行計(jì)算，該求解器所實(shí)現(xiàn)的動網(wǎng)格(ALE算法)技術(shù)可用來模擬推進(jìn)劑燃燒及變形帶來的拓?fù)浣Y(jié)構(gòu)改變的情況。在采用有限體積法進(jìn)行方程離散時(shí)，采用了高階的類WENO法，即HLLC法，該方法可處理類似點(diǎn)火流動等強(qiáng)瞬態(tài)過程。時(shí)間積分方法可通過3階或4階顯式多步Runge-Kutta時(shí)間步進(jìn)算法。

這兩個(gè)流體動力學(xué)求解器都可以選擇是否包含有湍流(Rocturb)、拉格朗日粒子追蹤(Rocpart)、歐拉煙霧追蹤(Rocsmoke)、化學(xué)反應(yīng)(Rocspecies)以及輻射(Rocrad)等模型。

2.4.2 結(jié)構(gòu)求解器

Rocstar主要有兩個(gè)有限元結(jié)構(gòu)力學(xué)求解器，Rocfrac和Rocsolid[13]。這兩個(gè)求解器中都包含有ALE算法，因此都可以考慮固體推進(jìn)劑由固相到氣相的變化過程。此外，還可以處理大的應(yīng)變、旋轉(zhuǎn)過程、三維熱傳導(dǎo)等問題。

Rocfrac是一個(gè)計(jì)算結(jié)構(gòu)力學(xué)和瞬態(tài)熱傳導(dǎo)求解器，該求解器采用非結(jié)構(gòu)網(wǎng)格，可進(jìn)行顯式和隱式求解。Rocfrac支持多種材料力學(xué)模型，主要為線彈性模型、Arruda-Boyce模型及Neoh ookian模型。Rocfrac在普通的四面體單元之間還可以包含有結(jié)合力/體有限單元，采用這種單元可以對藥柱裂紋的發(fā)展過程進(jìn)行仿真計(jì)算。

Rocsolid是一個(gè)有限應(yīng)變結(jié)構(gòu)求解器，它采用隱式時(shí)間積分方案，并且包含有多重網(wǎng)格和穩(wěn)定雙共軛梯度在內(nèi)的多個(gè)高性能稀疏迭代求解器，可以進(jìn)行高效并行線性系統(tǒng)的求解。對可壓縮或不可壓縮材料Rocsolid均包含多種求解單元類型。Rocsolid可以采用廣義有限元法(GFEM)來捕捉不連續(xù)性，例如運(yùn)動前緣或裂紋，而不需要繁瑣的幾何問題計(jì)算。除此之外，Rocsolid可利用多個(gè)基于微觀力學(xué)的模型來解釋微觀結(jié)構(gòu)的演變，包括裂紋形成和生長的影響，損傷(脫濕)以及應(yīng)變硬化。

2.4.3 燃燒求解器

Rocburn可仿真推進(jìn)劑加熱、點(diǎn)火及燃燒情況下燃面的退移速率或化學(xué)反應(yīng)。Rocburn有3個(gè)燃燒求解模型用于計(jì)算推進(jìn)劑燃率，物理模型在形式上是一維(沿燃面法向方向)，但是在推進(jìn)劑燃面上的每個(gè)單元面都獨(dú)立的施加，使其相當(dāng)于三維模型。這3個(gè)模型中最簡單的模型稱為RocburnAPN，采用著名的維也里公式r=apn，其燃速大小正比于當(dāng)?shù)厝細(xì)鈮簭?qiáng)的n次方。另外兩個(gè)燃速模型為動態(tài)燃速模型，通過求解一維非定常熱傳導(dǎo)方程以獲得溫度場從而捕捉點(diǎn)火瞬間過程。其中一個(gè)動態(tài)燃速模型稱為RocburnZN，基于Zeldovich-Novozhilov方法開發(fā)；另一個(gè)稱為RocburnPY，采用一個(gè)簡單的分解定律來計(jì)算燃速。RocburnPY可以計(jì)算點(diǎn)火器產(chǎn)生的高溫氣體對推進(jìn)劑的熱傳導(dǎo)過程以及推進(jìn)劑的點(diǎn)火過程。在一個(gè)全系統(tǒng)仿真里，RocburnPY可以調(diào)用由CSAR開發(fā)的詳細(xì)的三維推進(jìn)劑燃燒仿真代碼來計(jì)算推進(jìn)劑的燃速，此時(shí)可考慮推進(jìn)劑配方對燃速的影響。

Rocfire 用于模擬尺度在1～10 μm的異質(zhì)復(fù)合推進(jìn)劑的燃燒，異質(zhì)推進(jìn)劑顆粒填充幾何模型采用Rocpack生成。Rocfire考慮固體和氣相之間詳細(xì)的物理耦合，可考慮的現(xiàn)象包括燃面的非均勻退移、固體中的非穩(wěn)態(tài)熱傳導(dǎo)、金屬顆粒的存在、粘合劑/氧化劑的分解與火焰。通過Rocfire的模擬結(jié)果可提供氣相火焰結(jié)構(gòu)、燃速特性、渦度和湍流的產(chǎn)生的信息。

2.5 后處理

CSAR開發(fā)了一款功能強(qiáng)大的科學(xué)可視化后處理軟件包Rocketeer[29]。Rocketeer基于VTK開發(fā)，可在AIX、Linux、Windows及Solaris等系統(tǒng)上運(yùn)行。Rocketeer有一個(gè)非常友好的圖形界面，可用其讀取HDF文件格式或CGNS文件格式數(shù)據(jù)。Rocstar還開發(fā)了Rocketeer的批處理并行版本Voyager，以及客戶端/服務(wù)器版本Apollo/Houston。

除使用CSAR開發(fā)的專用后處理Rocketeer外，還可以使用商業(yè)軟件Tecplot和開源軟件ParaView對仿真結(jié)果進(jìn)行處理分析。

3 耦合方案

Rocstar的特點(diǎn)之一就是可處理三維情況下具有共同交接面的多域耦合問題。Rocstar要求各域可變形或移動但不能分開，也不能重疊。在交接面處，質(zhì)量、動量及能量等保持守恒，一些物理量，如密度等被設(shè)置為跳躍邊界條件。

Rocstar采用分域計(jì)算策略，不同計(jì)算域單獨(dú)進(jìn)行網(wǎng)格劃分，對于不同的物理求解域都有對應(yīng)的求解器，比如流體域的Rocflo和固體域的Rocfrac。每個(gè)域的計(jì)算過程是獨(dú)立的，只在每個(gè)全局時(shí)間步長計(jì)算后，將交接面的數(shù)據(jù)通過Rocface模塊傳遞給相鄰域求解器。例如流體域在計(jì)算一個(gè)全局時(shí)間步長后，將壓力分布加載到交接面上并傳遞給固體域，固體域進(jìn)行一個(gè)全局時(shí)間步長的計(jì)算，得到交接面的變形及運(yùn)動速度，又將其傳遞給流體域，整個(gè)過程一直重復(fù)直到計(jì)算結(jié)束為止，如圖5所示。交接面處可考慮燃燒或化學(xué)反應(yīng)過程，由燃燒或變形造成的計(jì)算域幾何形狀發(fā)生極端變化時(shí)，需要對網(wǎng)格進(jìn)行光順或重構(gòu)處理。

圖5 Rocstar多物理接口耦合Fig.5 Rocstar multiphysics interface coupled

Rocstar 采用“分割”方法來處理多域耦合時(shí)間步進(jìn)[23，30]。在這個(gè)方法中，一個(gè)全局時(shí)間步長Δtn內(nèi)，流體域的時(shí)間步長Δtf和固體域的時(shí)間步長Δts是不同的，由CFL條件來控制其大小。當(dāng)前域物理求解器在達(dá)到了下一全局時(shí)刻tn+1時(shí)，會通過交界面與相鄰域之間進(jìn)行數(shù)據(jù)傳遞。只要全局時(shí)間不比任何一個(gè)物理計(jì)算求解器的時(shí)間步長大太多，整個(gè)系統(tǒng)仍保持緊密的耦合。圖6為Rocstar的顯式時(shí)間步進(jìn)耦合方案。

圖6 Rocstar耦合時(shí)間步進(jìn)方案Fig.6 Rocstar coupled time stepping scheme

如圖6所示，在Δtn時(shí)刻計(jì)算首先從固體域開始，當(dāng)固體域求解器經(jīng)若干時(shí)間步Δts的計(jì)算到達(dá)下一全局時(shí)刻tn+1，將燃面的位置、速度以及由于燃燒帶來的質(zhì)量流量傳遞給流體域求解器，并采用Rocprop模塊更新燃面位置。流體域求解器經(jīng)過多個(gè)時(shí)間步Δtf的計(jì)算從tn時(shí)刻到達(dá)tn+1時(shí)刻。此時(shí)，流體域會將燃面處的受力加載傳給固體域，并將燃面處的壓強(qiáng)及溫度傳遞給燃燒求解器，燃燒求解器計(jì)算出的燃速值也會傳回給固體域求解器，這樣就完成了一個(gè)全局時(shí)間步的計(jì)算。

對于結(jié)構(gòu)求解，如果采用隱式固相結(jié)構(gòu)求解器Rocsolid代替Rocfrac，還可以采用基于預(yù)測-修正的隱式時(shí)間步算法。

4 Rocstar程序的驗(yàn)證

CSAR在開發(fā)Rocstar過程中，選擇航天飛機(jī)可重復(fù)使用固體助推器作為其主要仿真對象，NASA向CSAR提供該發(fā)動機(jī)的各種數(shù)據(jù)，以供CSAR的研究人員對Rocstar進(jìn)行各類驗(yàn)證[22]。

在開發(fā)Rocstar過程中，為了對各模塊及物理求解器進(jìn)行驗(yàn)證，CSAR開發(fā)人員進(jìn)行了相應(yīng)的數(shù)值實(shí)驗(yàn)對質(zhì)量守恒、網(wǎng)格精度、系統(tǒng)時(shí)間步和時(shí)間精度，以及ALE算法的可擴(kuò)展性進(jìn)行了研究，在開發(fā)物理求解器過程中還將計(jì)算結(jié)果與其他求解器進(jìn)行對比以保證Rocstar的計(jì)算準(zhǔn)確性與精度[31-33]。

除了航天飛機(jī)固體助推器發(fā)動機(jī)，CSAR還對多臺戰(zhàn)略火箭發(fā)動機(jī)進(jìn)行了仿真計(jì)算，證明了其仿真結(jié)果具有很高的可信度。例如，CSAR開發(fā)人員使用Rocstar對1991年美國“大力神4”助推固體火箭發(fā)動機(jī)(Titan IV SRMU)在首發(fā)試車時(shí)發(fā)生爆炸的原因進(jìn)行了數(shù)值模擬[34-35]。通過計(jì)算發(fā)現(xiàn)，該發(fā)動機(jī)是由于燃?xì)庋厮幹ǖ懒鲃哟嬖趬簭?qiáng)分布，導(dǎo)致后段藥柱變形，在連接槽處推進(jìn)劑發(fā)生塌陷，從而減少了燃?xì)饬鲃拥挠行饬髅娣e，加速了壓強(qiáng)分布的不均勻性，而這又加速導(dǎo)致藥柱的變形，如圖7所示，仿真結(jié)果與理論推斷相吻合[36]。

圖7 變形推進(jìn)劑中空隙的分布Fig.7 Distribution of voids in the deformed propellant

5 Rocstar仿真求解流程

Rocstar的仿真流程與多數(shù)CAE仿真軟件(如CFX、Fluent和OpenFOAM等)的仿真流程類似，分為前處理、運(yùn)算求解和后處理三部分。本節(jié)以Rocstar自帶的姿態(tài)控制發(fā)動機(jī)ACM(Attitude Control Motor)算例為例介紹采用Rocstar進(jìn)行固體火箭發(fā)動機(jī)仿真時(shí)，其仿真求解流程及相關(guān)模型、參數(shù)的設(shè)置方法。

該姿控發(fā)動機(jī)的構(gòu)型如圖8所示，其總長度約為55 mm，是一個(gè)小型的固體火箭發(fā)動機(jī)[33]。采用Rocstar對該發(fā)動機(jī)進(jìn)行非定常計(jì)算，在考慮燃面動態(tài)退移情況下，開展發(fā)動機(jī)內(nèi)流場燃燒流動耦合仿真計(jì)算。將發(fā)動機(jī)內(nèi)部流動簡化無粘流動，在推進(jìn)劑燃面采用質(zhì)量流量入口，推進(jìn)劑燃速采用燃速指數(shù)公式r=apn進(jìn)行求解。

圖8 姿態(tài)控制發(fā)動機(jī)構(gòu)型Fig.8 Attitude control motor geometry

5.1 前處理

采用幾何造型軟件Pro/Engineer繪制發(fā)動機(jī)幾何模型并導(dǎo)出通用幾何模型接口IGES格式文件。在Gridgen中，采用多塊結(jié)構(gòu)六面體網(wǎng)格對發(fā)動機(jī)進(jìn)行劃分(如圖9所示)，并輸出對應(yīng)的網(wǎng)格文件及邊界條件文件。

圖9 發(fā)動機(jī)多塊結(jié)構(gòu)網(wǎng)格Fig.9 Block-structured hexahedral mesh of motor

圖10 本地?cái)?shù)據(jù)存檔目錄結(jié)構(gòu)Fig.10 Set up NDA scheme

Rocstar沒有圖形操作界面，相關(guān)設(shè)置主要通過修改對應(yīng)的配置文件進(jìn)行。建立如圖10所示的NDA(Native Data Archive)目錄結(jié)構(gòu)[37]，其中每一個(gè)矩形表示一個(gè)文件夾，箭頭指向?yàn)楫?dāng)前文件夾的子文件夾。將相應(yīng)的網(wǎng)格文件(ACM-PLOT3D.grd)、邊界條件文件(ACM-PLOT3D.inp /ACM.bc)、邊界條件映射文件(ACM-PLOT3D.bcmp)、物理模型和數(shù)值參數(shù)等信息配置文件(ACM.inp)以及各求解模塊的參數(shù)文件(*.txt)放入相應(yīng)的位置。

當(dāng)上述工作準(zhǔn)備完成后采用2.3節(jié)中介紹的前處理工具Rocprep，為每個(gè)物理仿真應(yīng)用創(chuàng)建數(shù)據(jù)輸入集和進(jìn)行并行運(yùn)算分區(qū)劃分。最終生成的姿控發(fā)動機(jī)算例目錄結(jié)構(gòu)如圖11所示，其中ACM40為姿控發(fā)動機(jī)算例求解的主目錄。在ACM40目錄下包含3個(gè)子文件夾：Rocflo文件夾，即RocburnAPN文件夾和Rocman文件夾，用于存放對應(yīng)模塊的配置參數(shù)文件、相關(guān)數(shù)據(jù)信息文件和計(jì)算結(jié)果輸出文件等。

圖11 姿態(tài)控制發(fā)動機(jī)算例目錄結(jié)構(gòu)Fig.11 Directory structure of ACM case

5.2 運(yùn)算求解

前處理設(shè)置完成后，在運(yùn)算求解階段主要通過對圖11中各求解模塊的配置文件進(jìn)行參數(shù)設(shè)置。例如在RocstarControl.txt文件中，可以對耦合方案、流體求解模塊、燃燒求解模塊以及時(shí)間步長、自動保存步數(shù)等進(jìn)行設(shè)置。RocburnAPN文件夾內(nèi)的RocburnAPNControl.txt文件用于設(shè)置推進(jìn)劑燃速指數(shù)公式r=apn中的系數(shù)n和壓強(qiáng)指數(shù)a，以及絕熱火焰溫度等參數(shù)。

Rocstar采用MPI進(jìn)行大規(guī)模并行求解計(jì)算，在求解模型及參數(shù)設(shè)定完成后采用如下命令開始求解計(jì)算：mpirun-np核數(shù)rocstar。計(jì)算結(jié)果保存在對應(yīng)模塊文件夾內(nèi)的Rocout文件夾。

5.3 后處理

計(jì)算完成后，采用Rocketeer后處理軟件對仿真結(jié)果進(jìn)行處理分析。圖12為選取計(jì)算結(jié)果中4個(gè)時(shí)刻下發(fā)動機(jī)幾何和網(wǎng)格拓?fù)渥兓?，可看出隨著推進(jìn)劑的燃燒，燃面發(fā)生了退移。圖13為z=0平面發(fā)動機(jī)內(nèi)流場壓強(qiáng)云圖，可以看出隨著推進(jìn)劑的燃燒，發(fā)動機(jī)內(nèi)流道變大且發(fā)動機(jī)內(nèi)部壓強(qiáng)逐漸升高。

圖12 幾何構(gòu)型及網(wǎng)格拓?fù)渥兓疐ig.12 Topology change of geometry and mesh

圖13 Z=0平面壓強(qiáng)變化Fig.13 Pressure change at Z=0 plane

6 結(jié)束語

Rocstar多物理場仿真程序提供了先進(jìn)的固體火箭發(fā)動機(jī)仿真模擬能力，在固體火箭發(fā)動機(jī)全三維模型和精細(xì)物理現(xiàn)象捕捉、多物理場非線性耦合、復(fù)雜幾何模型動態(tài)變化、空間和時(shí)間尺度上的級大差異、材料屬性以及物理過程的復(fù)雜性以及大規(guī)模并行實(shí)現(xiàn)等方面都實(shí)現(xiàn)了突破。Rocstar能夠?qū)腆w火箭發(fā)動機(jī)工作過程進(jìn)行全三維、多域耦合、高精度且大規(guī)模并行的數(shù)值模擬，可用于預(yù)示固體火箭發(fā)動機(jī)的性能和工作效果，分析潛在可能導(dǎo)致發(fā)動機(jī)工作失效的問題，以及為事故原因分析和發(fā)動機(jī)結(jié)構(gòu)設(shè)計(jì)的合理改進(jìn)提供較準(zhǔn)確的技術(shù)支持。

Rocstar程序目前已經(jīng)開源，除采用Rocstar對固體火箭發(fā)動機(jī)開展已有模型算法的仿真計(jì)算，還可以對其源程序代碼進(jìn)行二次開發(fā)。例如，Rocstar中的兩相流中并沒有考慮顆粒之間的碰撞過程，如需在仿真過程中添加該過程，則可以對Rocpart拉格朗日顆粒追蹤模塊進(jìn)行修改。此外，得益于其靈活、通用的集成框架和眾多的求解模塊，還可以采用Rocstar對固體火箭沖壓發(fā)動機(jī)或固體火箭超燃沖壓發(fā)動機(jī)的工作性能開展仿真計(jì)算。