自適應間斷裝配法模擬彈道靶中超高速彈丸發(fā)射

常思源, 鄒東陽, 劉 君,*

(1. 大連理工大學 航空航天學院, 遼寧 大連 116024; 2. 中國空氣動力研究與發(fā)展中心, 四川 綿陽 621000)

0 引 言

超聲速自由飛彈道靶是專門用于炮彈、火箭、導彈或其他飛行器模型的自由飛試驗、進行氣動力和彈道性能研究的封閉式射擊靶道[1]，是現(xiàn)代氣動試驗最重要的設(shè)備之一。模型發(fā)射裝置是彈道靶的重要組成部分。當前使用最廣泛的模型發(fā)射裝置是二級輕氣炮，其工作原理是利用火藥爆炸對泵室內(nèi)的輕氣(如氫氣或氦氣)進行加溫加壓，驅(qū)動模型(彈丸與彈托)在較長的發(fā)射管中加速到極高的速度(如5km/s)[2]。

利用彈道靶發(fā)射彈丸時，必須借助彈托分離技術(shù)[3]。通常將彈丸固定在用塑料制成的彈托內(nèi)，彈托外徑與發(fā)射管內(nèi)徑相同；加速過程中，彈托與發(fā)射管接觸摩擦；出膛后，利用氣動力分離技術(shù)使彈托逐漸偏離預定彈道并被攔截器阻擋，僅使彈丸進入試驗段開展測試。

當前，國內(nèi)面向彈道靶開展的數(shù)值研究主要集中于運用經(jīng)典內(nèi)彈道理論以及計算流體力學(CFD)方法對二級輕氣炮開展相關(guān)研究，為優(yōu)化二級輕氣炮性能提供可靠的數(shù)值分析手段[4-6]。彈丸從靜止啟動到加速、出膛、分離的過程，涉及的物理過程復雜、空間尺度差異較大、非定常效應明顯、激波強度較高，對數(shù)值模擬提出了一定挑戰(zhàn)，調(diào)研發(fā)現(xiàn)相關(guān)研究國內(nèi)鮮有報道。

彈道靶中彈丸發(fā)射全過程存在強烈的氣動干擾，強激波、膛口波系、高速膨脹氣體都會對彈丸飛行運動造成比較大的影響。發(fā)展彈丸發(fā)射全過程數(shù)值模擬技術(shù)有助于彈道靶試驗設(shè)計，以保證彈丸運動軌跡的準確性，避免出現(xiàn)彈丸偏離預定彈道而被攔截或打到靶室洞壁的情況。

本文采用劉君[7-8]等基于非結(jié)構(gòu)動網(wǎng)格技術(shù)和格心型有限體積方法發(fā)展的自適應間斷裝配求解器(Adaptive Discontinuity Fitting solver，ADFs)對非定常流場中的運動激波進行裝配，通過二維算例對彈丸發(fā)射全過程開展數(shù)值仿真，較為詳細地介紹各個運動階段數(shù)值模擬的特點。

1 計算方法

1.1 控制方程

計算主要考慮激波的影響，未考慮粘性效應?？刂品匠滩捎没贏LE(Arbitrary Lagrange-Euler)方法的可壓縮非定常Euler方程的積分形式：

(1)

式中，Q為守恒變量，F(xiàn)c為對流項通量，xc為網(wǎng)格運動速度，Ω、?Ω分別為控制體和控制邊界，n為控制體邊界外法向向量，dV、dS分別為體積微元和面積微元。

使用格心型有限體積方法時，流動變量存儲在計算網(wǎng)格單元的中心。根據(jù)單元界面兩側(cè)的參數(shù)，通過通量矢量分裂方法計算對流項通量Fc，如圖1所示。若使用左右單元的格心值代表界面兩側(cè)的參數(shù)，則空間精度僅有一階。為獲得二階精度，使用Gauss-Green公式計算單元梯度[9]：

(2)

(3)

其中，Np為第k個面元所包含的節(jié)點數(shù)，Nc是與第n個網(wǎng)格節(jié)點相鄰的單元的總數(shù)，cm是第m個鄰居單元的距離反比權(quán)系數(shù)，即cm=1/|rmn|。

圖1 流場變量重構(gòu)

1.2 非結(jié)構(gòu)動網(wǎng)格技術(shù)

高速條件下的多體分離流場呈現(xiàn)高度非線性和非定常的特征。對此類問題進行數(shù)值模擬，按照依賴網(wǎng)格技術(shù)的差異主要分為3種方法，即動態(tài)嵌套網(wǎng)格[10-11]、非結(jié)構(gòu)動網(wǎng)格[12-13]和笛卡爾自適應網(wǎng)格[14-15]。本文使用基于頂點模型彈簧近似法的非結(jié)構(gòu)動網(wǎng)格技術(shù)[16]，下面對該技術(shù)進行簡要介紹。

如圖2所示，彈簧近似法將網(wǎng)格單元的每條邊視作彈簧，認為網(wǎng)格點i和j之間的彈簧張力與兩節(jié)點之間距離成正比，滿足Hooke定律。因此，對作用于節(jié)點i的所有彈性力進行合成，可得：

(4)

其中，Ni是與節(jié)點i相連的節(jié)點總數(shù)，Kij是網(wǎng)格邊ij的彈簧倔強系數(shù)。

圖2 彈簧近似模型示意圖

根據(jù)初始時刻的網(wǎng)格，按照上式可以得到初始狀態(tài)下每個網(wǎng)格點i受到的合力。當某些節(jié)點移動后，彈簧長度變化使得該點的受力發(fā)生變化，為了保持整體受力平衡，需要調(diào)整其他網(wǎng)格點的位置。由于網(wǎng)格變形過程中每點的合力不變，新的網(wǎng)格分布依然滿足上式。

1.3 激波裝配技術(shù)

對于含有強激波的流動問題，若單純使用捕捉方法進行處理，會存在數(shù)值不穩(wěn)定問題，同時也會影響流場的計算精度[17]。計算激波一個最直接的辦法就是使用激波裝配方法，即將激波面看作一個間斷，使兩側(cè)的物理量在激波法向上滿足Rankine-Hugoniot(R-H)公式：

(5)

式中，ρ1和ρ2是激波上下游密度，u1和u2是激波上下游的流動速度，p1和p2是上下游壓強，h1和h2是上下游比焓。

最早的激波裝配方法可以追溯到20世紀40年代，經(jīng)過半個多世紀的研究，特別是近十幾年的發(fā)展，其實用性已有了明顯提高[18]。劉君等人基于非結(jié)構(gòu)動網(wǎng)格技術(shù)發(fā)展出的ADFs具有較強靈活性，能夠很好地滿足本文模擬非定常流動下強激波的需求。

如圖3所示，在二維情況下，對于初始流場中預先給定的激波線，其上各個網(wǎng)格節(jié)點的屬性可以在激波節(jié)點(S)和普通節(jié)點(O)之間靈活轉(zhuǎn)換，自動適配當前流場狀態(tài)；即，當實際激波尚未到達給定線時，該線不作為激波進行裝配計算，僅當實際激波到達時，才轉(zhuǎn)化為激波進行裝配。在下文2.1節(jié)有相關(guān)具體闡述。

圖3 網(wǎng)格節(jié)點屬性定義

相對于定常計算，在非定常計算過程中，運動激波的裝配無需進行特殊處理，這也是ADFs的一個優(yōu)點。具體來講，當推進到n+1時刻時，首先根據(jù)n時刻物理量通過R-H關(guān)系式迭代求解各個激波節(jié)點的運動速度ω(如圖4所示，對于定常計算收斂時，ω→ 0)，從而得到n+1時刻激波陣面位置；隨后使用前文介紹的網(wǎng)格變形算法更新內(nèi)部網(wǎng)格節(jié)點的位置，當網(wǎng)格變形較大時，需要進行網(wǎng)格重構(gòu)；最后進行時間推進計算n+1時刻的物理量。也就是說，在原有的激波捕捉法求解器基礎(chǔ)上，ADFs僅增加了一個裝配激波的模塊，通過輸出激波節(jié)點參數(shù)(物理量及激波速度)和激波面元通量，建立了與捕捉法求解器的聯(lián)系，故該算法很容易移植到任何含有動網(wǎng)格的程序中。限于篇幅，對ADFs的詳細闡述請參閱文獻[19]。

圖4 激波節(jié)點參數(shù)確定示意圖

在空間離散上，采用目前流行的van Leer矢通量分裂格式，使用前文介紹的梯度重構(gòu)方法，保證了空間離散具有二階精度；時間推進采用顯式四階四步Runge-Kutta方法，時間步長的量級為10-6；使用商業(yè)軟件POINTWISE劃分全場網(wǎng)格，采用Advancing Front算法保證絕大部分區(qū)域網(wǎng)格為正三角形，提高了整體網(wǎng)格的質(zhì)量。

2 數(shù)值仿真

在一些彈道靶試驗中，通過點燃火藥后形成的高溫高壓燃氣驅(qū)動活塞對泵室內(nèi)的氣體加壓；壓力達到一定數(shù)值后，泵室和發(fā)射管之間的膜片破裂，高壓氣體直接作用于彈托底部，驅(qū)動彈丸/彈托一起加速運動；彈丸/彈托在比較長的發(fā)射管道中加速至預定速度，進入分離罐，在發(fā)射管出口的高壓氣體和彈丸/彈托附近的激波共同作用下，經(jīng)過合理設(shè)計的彈托和彈丸可以實現(xiàn)氣動力分離。

計算模型如圖5所示，主要包括發(fā)射管和分離罐兩部分。發(fā)射管總長10m，管道內(nèi)徑0.06m；彈丸呈子彈頭形，長0.10m，底部直徑0.03m。彈丸在兩個完全相同的彈托約束下沿著發(fā)射管水平向右方運動。根據(jù)文獻[5]的數(shù)據(jù)，直接給定彈丸和彈托(后文也將彈丸和彈托合稱為“模型”)在發(fā)射管內(nèi)的加速過程，其沿x方向(發(fā)射管軸向)的運動速度vx和位移L的曲線如圖6所示。

圖5 彈道靶計算模型示意圖

Fig.5Schematicdiagramofthecomputationalmodelforballisticrange

根據(jù)該問題的特點，整個試驗過程可分為模型加速、激波出膛、模型出膛、氣動分離4個階段。

圖6 模型運動速度與位移曲線

2.1 模型加速過程

該過程的一個顯著特點是：物體的空間尺度相差較大，即與加速距離10m相比，模型尺寸要小若干個量級。如果直接對整體發(fā)射管進行建模，采用激波捕捉法求解模型沿發(fā)射管的運動，就會付出相當大的計算代價，對于三維問題而言尤其難以接受。

分析發(fā)現(xiàn)，模型加速運動會壓縮處于上游的氣體，使其頭部某處形成一道向右運動的正激波，且激波陣面為模型加速運動形成的擾動的最遠傳遞位置，而激波陣面上游流體則始終處于未擾動的靜止狀態(tài)。由于模型加速過程給定，可以不考慮其底部流動狀態(tài)的變化，即在該過程中僅需計算激波和模型之間區(qū)域的流動變化。

圖7給出了初始計算區(qū)域(x方向長約0.1m，占整個發(fā)射管長度的1/100；Δx為彈丸頭部頂點與運動激波的水平距離)，以非結(jié)構(gòu)三角形網(wǎng)格進行空間離散。設(shè)置模型右側(cè)邊界為激波邊界，當激波未形成或激波未達到該位置時，使用ADFs可以切換激波邊界為出口邊界，不進行裝配，此時計算網(wǎng)格隨模型一起剛性運動，不發(fā)生變形；當激波形成且到達邊界位置時，開始采用裝配方法進行處理，激波邊界條件啟動。

圖7 初始計算區(qū)域和邊界條件

圖8給出了加速階段5個時刻的壓力云圖?？梢钥闯?，由于激波運動速度大于模型運動速度，計算過程中計算區(qū)域會逐漸擴大。圖9給出了3個時刻激波附近的網(wǎng)格分布，自激波下游某處(圖9中四邊形網(wǎng)格和三角形網(wǎng)格交界處)起至模型的這部分網(wǎng)格，可以設(shè)置為隨模型一起運動的剛性運動網(wǎng)格，使全場網(wǎng)格變形及網(wǎng)格重構(gòu)區(qū)僅包含激波附近若干層網(wǎng)格，從而提高全場網(wǎng)格質(zhì)量，減少網(wǎng)格重構(gòu)次數(shù)，極大地提高計算效率。

圖8 加速階段不同時刻的壓強云圖

Fig.8Pressureiso-contoursatdifferentmomentsintheacceleratedprocess

圖9 不同時刻激波附近的網(wǎng)格分布(從左至右：t=2.00ms, 3.00ms, 4.00ms)

Fig.9Meshdistributionsneartheshockatdifferentmoments(fromlefttoright:t=2.00ms,3.00ms,4.00ms)

圖10為模型及激波的加速過程運動特性圖(圖中Mas為運動激波的馬赫數(shù))，圖11為彈丸和上彈托迎風面在加速階段所受氣動力F的變化曲線(Fx、Fy分別表示x、y方向的氣動力)。在初始加速階段，激波相對模型的運動速度較小，但加速度較高；在t=1.00ms時刻左右，激波運動速度超過模型運動速度；隨后兩者的速度差逐漸增大，激波強度逐漸增強，模型表面氣動力開始快速升高；當模型運動加速度發(fā)生改變后(t=2.00ms)，激波的運動加速度也隨之改變。因此，如果該過程的加速度曲線不同，即使模型到達發(fā)射管出口的速度一樣，不同時刻的激波強度以及模型表面的壓強也會有所差異。

圖10 模型與激波的運動特性曲線

圖11 加速過程彈丸和彈托的受力曲線

Fig.11Aerodynamicforcecurvesofprojectileandsabotsintheacceleratedprocess

值得注意的是：在模型加速過程中，激波前的相對馬赫數(shù)會超過15，若在非結(jié)構(gòu)網(wǎng)格中使用常規(guī)捕捉方法計算，很容易出現(xiàn)數(shù)值震蕩，對計算方法的穩(wěn)定性造成影響；而使用裝配方法則能避免這個問題，這也是使用裝配方法進行計算的另一個優(yōu)點。

2.2 激波出膛過程

激波出膛過程從右行運動激波到達發(fā)射管口即將進入分離罐時刻開始計算，直到模型即將出膛終止計算。此階段激波在分離罐中逐漸向前推進，最終形成了一個不斷外擴的右行弧狀激波陣面。此過程仍使用激波裝配法處理該激波，計算區(qū)域隨時間逐漸擴大。

圖12給出了5個時刻激波陣面和模型的位置與形狀變化，圖13為中間某時刻整個流場的壓強與馬赫數(shù)分布云圖。計算初始時刻(t=4.47ms)，在發(fā)射管出口預設(shè)一個擬裝配激波面，當激波未到達該位置時，流動按照捕捉方法進行計算；激波到達后，該界面逐漸全部變?yōu)榧げ嚸孢M行裝配計算。

激波出膛之后，在分離罐壁面的作用下膨脹形成一個半圓形的影響區(qū)域。由于氣流膨脹引起下游壓強降低，靠近發(fā)射管軸線處的激波運動速度迅速降低，最終會小于模型運動速度；而在靠近壁面處，由于氣流過膨脹，壓力降幅更大，激波運動速度更小，最終造成靠近壁面的激波陣面向下游偏斜。

圖12 不同時刻激波陣面的位置

圖13 t=4.55ms時刻流場壓強與馬赫數(shù)云圖

2.3 模型出膛過程

模型出膛過程涉及到動態(tài)多體分離中一個比較常見的問題——“接觸-分離問題”。

在發(fā)射出膛過程中，彈丸/彈托與彈托/發(fā)射管之間都保持接觸狀態(tài)，是一個典型的接觸-分離問題。該過程存在網(wǎng)格拓撲變化，很難使用一套網(wǎng)格進行描述。處理此類問題，往往需要人為在物體之間設(shè)置一定縫隙。為盡可能減小縫隙對流動造成的影響，文獻[20]發(fā)展出一種“虛擬網(wǎng)格通氣”技術(shù)(如圖14所示)，在高壓氣體和縫隙的網(wǎng)格交界面設(shè)置“虛擬擋板”，防止彈丸/彈托上下游高壓氣體通過縫隙相互影響，物理上實現(xiàn)物體間的緊密接觸狀態(tài)。

模型出膛后，其底部高壓氣體溢出，從而影響靶室內(nèi)的流動狀態(tài)。由于模型運動速度較快，相對時間尺度較小，為簡化處理，假設(shè)發(fā)射管內(nèi)高壓氣體的狀態(tài)不變，即p=60MPa,ρ=40kg/m3。

2.4 氣動分離過程

待完全飛出發(fā)射管，模型開始在氣動力作用下進行六自由度運動。耦合求解流動方程和六自由度運動方程，獲得流場以及模型位移和姿態(tài)的變化情況。在此過程中，可以使用前文介紹的非結(jié)構(gòu)動網(wǎng)格技術(shù)描述模型運動軌跡。

圖15為分離過程中5個時刻的馬赫數(shù)云圖，圖16、17分別為彈丸、彈托的速度vx以及氣動力F的變化曲線。從云圖中可以看出，在模型上游兩道激波和尾部高速高壓氣流共同建立的流場下，模型逐漸追上前方兩道激波。彈丸、彈托分別在t=4.84ms和4.87ms時刻開始穿過第一道內(nèi)部膛口激波，進入激波下游高壓區(qū)后，氣動力急劇增大，x方向的運動速度逐漸降低。同時，由于彈托外形的不對稱設(shè)計，其在y方向受到很大的氣動力，迅速偏離彈丸的運動軌跡。在t=4.95ms時刻，彈丸即將追上并穿過初始右行弧狀激波，從而擺脫該激波流場，進入試驗段開展相關(guān)測試。

圖14 模型出膛過程的“虛擬擋板”示意圖

Fig14Schematicofvirtualbafflesintheprocessofthemodelleavingthemuzzle

圖15 分離階段不同時刻馬赫數(shù)云圖

Fig.15Machiso-contoursatdifferentmomentsintheprocessofsabotseparation

圖16 彈丸和彈托速度變化曲線

圖17 彈丸和彈托氣動力變化曲線

需要說明的是，雖然對發(fā)射過程開展二維模擬與實際情況不相符合，但是所得到的氣動力可以定性地對這一過程的流動變化進行說明。

3 結(jié) 論

將彈道靶試驗中以超高速發(fā)射的彈丸作為研究對象，通過對二維問題的模擬，系統(tǒng)介紹了數(shù)值模擬時所需要的一些技術(shù)手段，如剛性和變形網(wǎng)格技術(shù)、虛擬網(wǎng)格通氣技術(shù)、激波裝配方法等。一方面，采用基于非結(jié)構(gòu)動網(wǎng)格技術(shù)和格心型有限體積方法發(fā)展的自適應間斷裝配法，對該流場中大范圍運動的高強度激波進行了裝配計算，拓展了激波裝配方法在工程問題中的應用；另一方面，針對彈道靶中超高速彈丸發(fā)射這類問題，建立了一套高效的數(shù)值模擬方法，實現(xiàn)了對彈丸從靜止啟動到加速、出膛、分離的全動態(tài)過程的精細流場刻畫，有助于指導試驗方案驗證和優(yōu)化設(shè)計。

一般而言，驅(qū)動段高壓氣體(如氫氣)從發(fā)射管溢出后與空氣相遇，還可能存在激波誘導燃燒等問題，因此，在下階段工作中需要進一步考慮流動中的化學反應等。