含化學(xué)反應(yīng)膛口流場的無網(wǎng)格數(shù)值模擬*

吳 偉,許厚謙,王 亮,薛 銳

(南京理工大學(xué)能源與動力工程學(xué)院，江蘇 南京 210094)

?

含化學(xué)反應(yīng)膛口流場的無網(wǎng)格數(shù)值模擬*

吳 偉,許厚謙,王 亮,薛 銳

(南京理工大學(xué)能源與動力工程學(xué)院，江蘇 南京 210094)

基于無網(wǎng)格方法，對包含大位移運動邊界和非平衡化學(xué)反應(yīng)的膛口流場進行了數(shù)值模擬。所發(fā)展算法是基于線性基函數(shù)最小二乘顯式無網(wǎng)格方法，忽略黏性及湍流的影響，對流場采用ALE(arbitrary Lagrangian-Eulerian)形式的Euler方程描述，對流通量和化學(xué)反應(yīng)源項采用多組分HLLC(Harten-Lax-van Leer-Contact)格式和有限速率反應(yīng)模型計算，對于運動邊界造成的點云畸形采用局部點云重構(gòu)方法處理，重構(gòu)過程中采用虛擬邊陣面推進。對圓柱繞流和激波誘導(dǎo)燃燒流場進行了數(shù)值模擬，驗證了重構(gòu)方法和化學(xué)反應(yīng)計算的有效性。最后對12.7 mm口徑機槍膛口流場進行了模擬，結(jié)果同實驗照片、非結(jié)構(gòu)網(wǎng)格方法結(jié)果吻合較好，數(shù)值結(jié)果清晰地再現(xiàn)了膛口初始沖擊波、膛口沖擊波、欠膨脹射流波系結(jié)構(gòu)的動力學(xué)發(fā)展過程，以及膛口焰的時間、空間分布特征。

爆炸力學(xué)；動態(tài)點云；無網(wǎng)格方法；膛口流場；非平衡反應(yīng)流

彈丸從膛口射出后，膛內(nèi)高溫、高壓的火藥氣體因突然被釋放而在膛口外急劇膨脹，形成氣動結(jié)構(gòu)異常復(fù)雜的膛口射流流場；負(fù)氧平衡的火藥氣體還可能與周圍的空氣發(fā)生劇烈的非平衡化學(xué)反應(yīng)，形成膛口焰；此外，流場幾何結(jié)構(gòu)也較為復(fù)雜并且包含高速運動的彈丸，這些都給膛口流場的數(shù)值模擬帶來了巨大的挑戰(zhàn)。目前，其數(shù)值研究大都建立在網(wǎng)格離散的基礎(chǔ)上，并且進行了不同程度的簡化，如簡化彈丸的形狀，忽略化學(xué)反應(yīng)、黏性等。

近幾十年，由于實際工程問題日益復(fù)雜，傳統(tǒng)網(wǎng)格方法出現(xiàn)了一些弊端，而新興的無網(wǎng)格方法得到了計算流體力學(xué)領(lǐng)域內(nèi)大量學(xué)者的關(guān)注，該方法擺脫了對網(wǎng)格的依賴，采用一系列節(jié)點離散求解域，通過構(gòu)建各節(jié)點的點云，采用最小二乘等方法直接求解強形式的控制方程。由于其求解過程只依賴點與點的聯(lián)系，無需構(gòu)造網(wǎng)格，因而對于復(fù)雜外形具有更強的適應(yīng)性，流場內(nèi)部布點也更為方便快捷。此外，相對于傳統(tǒng)網(wǎng)格方法，無網(wǎng)格方法最大的優(yōu)勢是易于處理流場中的大位移運動邊界，處理過程中均只需更新相應(yīng)點的點云信息即可，較非結(jié)構(gòu)網(wǎng)格效率更高。正是由于其內(nèi)在的靈活性、優(yōu)越性，無網(wǎng)格方法已經(jīng)得到了大量研究，并取得了一定成果。AUSM+-UP[1-2]、CUSP[3]、HLLC[4]等高精度、高激波分辨率格式成功應(yīng)用于無網(wǎng)格方法，對包含大位移運動邊界問題發(fā)展了重疊點云[5]和局部點云重構(gòu)[4]技術(shù)，此外還發(fā)展了一些無網(wǎng)格-笛卡爾網(wǎng)格混合網(wǎng)格算法[1,6-7]。近來，采用無網(wǎng)格方法對激波誘導(dǎo)燃燒流場進行模擬的算法[8]也見報道，但對包含大位移運動剛體的非平衡化學(xué)反應(yīng)流場還無法進行有效模擬。

本文中，采用無網(wǎng)格方法求解ALE形式多組分Euler控制方程，用于非平衡化學(xué)反應(yīng)流場的數(shù)值模擬，并結(jié)合局部點云重構(gòu)方法處理流場中的大位移運動邊界。首先，介紹改進的重構(gòu)方法的基本思想，以及所采用的無網(wǎng)格方法；隨后，對圓柱繞流和激波誘導(dǎo)燃燒流場進行模擬，驗證重構(gòu)方法和化學(xué)反應(yīng)計算的正確性；最后，對12.7 mm口徑機槍膛口流場進行模擬。

1 控制方程

對于包含任意位移運動邊界的非平衡化學(xué)反應(yīng)流場，忽略黏性及湍流的影響，采用ALE形式的多組分Euler方程描述：

(1)

式中：U為守恒變量，F(xiàn)、G為對流通量，W為化學(xué)反應(yīng)源項，S為軸對稱源項。

U、F、G、W具體定義如下：

2 點云重構(gòu)

對于運動邊界造成的點云畸形采用重構(gòu)的思想處理。首先,將相鄰的衛(wèi)星點以及各衛(wèi)星點與中心點互聯(lián)，形成虛擬邊；計算動邊界附近離散點點云質(zhì)量，查找質(zhì)量不滿足計算要求的離散點，進而形成點云重構(gòu)的空腔，避免了復(fù)雜的點云相交性判斷[4]；空腔內(nèi)采用虛擬邊推進的方法布點，同時生成網(wǎng)格信息(后處理軟件需要網(wǎng)格拓?fù)湫畔?，同填充布點[4]相比，該過程避免了后續(xù)點云Delaunay三角化過程，提高了效率；布點結(jié)束后更新空腔邊界離散點和新生成離散點的點云，并進行Laplace光順處理，重新計算其形函數(shù)；最后,采用線性插值的方法計算新生成離散點的物理量。

3 數(shù)值方法

3.1 空間離散

本文中發(fā)展的無網(wǎng)格方法采用線性基函數(shù)最小二乘擬合[4,9]直接求解各離散點的空間導(dǎo)數(shù)。假設(shè)流動基本變量滿足如下線性關(guān)系：

(2)

以任意點i(假設(shè)周圍分布6個衛(wèi)星點)為例，點i及其衛(wèi)星點均滿足式(2)，將中心點與各衛(wèi)星點進行疊加后可得：

(3)

令上述矛盾方程組的系數(shù)矩陣為A，采用最小二乘求解可得：

(4)

進而可得中心點i的空間導(dǎo)數(shù)為：

(5)

形函數(shù)bij、cij為(ATA)-1AT矩陣的第1、2行，進而中心點i的對流通量則可以表示為：

(6)

對中心點i與其衛(wèi)星點j的中點的通量Fij采用多組分HLLC格式[9]計算，并采用MUSCL重構(gòu)方法[1,4]對中點左右兩側(cè)狀態(tài)進行重構(gòu)，進而提高計算精度。應(yīng)用于無網(wǎng)格方法的HLLC格式具體形式如下：

(7)

p*=ρi(uni-Si)(uni-SM)+pi

(9)

un為流體速度的法向分量：

(10)

式中：v為流體的速度矢量。

波速Si、Sj以及接觸間斷面速度SM計算如下：

(11)

(12)

(13)

(14)

3.2 化學(xué)反應(yīng)模型

化學(xué)反應(yīng)源項采用有限速率反應(yīng)模型，反應(yīng)體系中的任意反應(yīng)可表示為：

(15)

(16)

式中：Am為指前因子，bm為溫度因子，Em為活化能，T為混合氣體溫度，Ru為通用氣體常數(shù)。逆向反應(yīng)速率常數(shù)Kbm可由反應(yīng)平衡常數(shù)計算。各離散點任意組分i的質(zhì)量生成率可由下式計算：

(17)

式中：NR為化學(xué)反應(yīng)總數(shù)，Mi為組分i的摩爾質(zhì)量。

3.3 時間項及邊界條件

采用四階Runge-Kutta法進行時間顯式推進。邊界采用的是在流場外構(gòu)造鏡像點的方法處理，鏡像點流動變量的取值則根據(jù)邊界類型確定。固壁采用法向無穿透邊界條件，遠場采用基于Riemann不變量的無反射邊界條件。

4 驗證算例

為驗證本文算法正確性，分別對圓柱繞流以及H2/Air預(yù)混氣體中的激波誘導(dǎo)燃燒流場進行了模擬，并與相關(guān)結(jié)果進行了比較。

4.1 圓柱繞流流場

為了驗證點云重構(gòu)技術(shù)的有效性，對圓柱繞流問題進行了模擬，計算域為2.2 m×1.4 m，圓柱半徑為0.1 m，初始時刻位于(0,0)處，計算分2種情形:(1)來流速度v0=0，圓柱運動速度vc=-600 m/s；(2)v0=600 m/s，vc=0。氧氣、氮氣摩爾比為1∶4，壓力p0、溫度T0分別為101.325 kPa、286.15 K，計算中凍結(jié)化學(xué)反應(yīng)。t=0.8 ms時的壓力云圖如圖1所示，圓柱上表面的量綱一壓力p/p0分布如圖2所示，可見采用動態(tài)點云重構(gòu)計算得到的結(jié)果同圓柱靜止時結(jié)果吻合較好。

圖1 壓力云圖(t=0.8 ms)Fig.1 Pressure contours at t=0.8 ms

圖2 圓柱上表面量綱一壓力分布Fig.2 Distribution of the dimensionless pressure along the upper surface of the cylinder

4.2 激波誘導(dǎo)燃燒流場

為驗證采用無網(wǎng)格方法計算非平衡化學(xué)反應(yīng)流的有效性，對激波誘導(dǎo)燃燒流場進行了模擬，該算例已經(jīng)成為驗證非平衡化學(xué)反應(yīng)流算法正確性的標(biāo)準(zhǔn)算例。自由來流為等化學(xué)當(dāng)量比的氫氣/空氣預(yù)混氣體，來流速度v0、壓力p0、溫度T0分別為1 685 m/s、42.662 kPa、250 K，化學(xué)反應(yīng)采用7組分8步反應(yīng)機理[11]。圖3為等溫線分布圖，可見本文算法成功捕捉到了流場中燃燒波與激波的分離，激波位置同實驗照片中的結(jié)果亦吻合較好。駐點流線上溫度、壓力以及主要組分質(zhì)量分?jǐn)?shù)的分布如圖4、5所示，同M.Soetrisno等[11]的結(jié)果基本一致。

圖3 溫度云圖Fig.3 Temperature contours

圖4 駐點流線上溫度、壓力分布Fig.4 Distribution of temperature and pressure along the stagnation stream line

圖5 駐點流線上組分質(zhì)量分?jǐn)?shù)Fig.5 Distribution of mass fractions along the stagnation stream line

圖6 膛口離散點分布示意圖 Fig.6 Schematic diagram for the distribution of the discrete points at the muzzle

5 膛口二次燃燒流場

對口徑為12.7 mm的高射機槍膛口流場進行了數(shù)值模擬，膛管的內(nèi)徑為13.7 mm，外徑為31.0 mm，膛管長為1.08 m，計算中對彈丸結(jié)構(gòu)進行了適當(dāng)?shù)暮喕?，其出膛速度?10 m/s。外流場取0.32 m×0.8 m矩形區(qū)域，流場共布點228 247個，某時刻流場局部區(qū)域離散點分布見圖6。定義彈底到達膛口瞬間為0時刻。僅在彈丸出膛后開啟化學(xué)反應(yīng)計算，其化學(xué)反應(yīng)采用H2-CO-O29組分11步反應(yīng)機理[10]。

彈丸在膛內(nèi)運動部分時刻的溫度云圖見圖7，在火藥氣體的作用下，彈丸由膛底開始運動，壓縮彈前空氣，形成一道具有一定強度的激波，彈丸與壓縮的空氣一同向膛口運動，到達膛口后向周圍靜止大氣環(huán)境迅速膨脹，形成包含初始沖擊波、馬赫盤、瓶狀激波、渦環(huán)等復(fù)雜氣動結(jié)構(gòu)的初始射流流場。

當(dāng)彈丸運動出膛口瞬間(t=0)，需對膛內(nèi)火藥氣體組分質(zhì)量密度、壓力、速度等物理量重新賦值，具體如下：

px=pd[1+φ(1-x2/L2)]

vx=v0x/L

(18)

式中：px、vx為膛管內(nèi)壓力和速度，x為距膛底距離，pd為彈丸運動至膛口時彈底壓力，取值74 MPa，φ取值0.18。火藥氣體平均密度ρ0為120 kg/m3，各組分的摩爾分?jǐn)?shù)分別為：CO，0.434 69；CO2，0.099 43；H2，0.135 81；H2O，0.222 63；N2，0.107 44。對溫度、總能根據(jù)熱力學(xué)關(guān)系計算獲得。

圖8為不同時刻計算陰影圖與實驗陰影照片的對比，可見t=0時刻初始沖擊波、t=250 μs時刻膛口沖擊波、馬赫盤的位置、形狀等同實驗結(jié)果吻合較好。圖9為采用本文算法計算得到的t=415 μs時刻密度、溫度云圖同非結(jié)構(gòu)網(wǎng)格方法[10]得到結(jié)果的對比。由于本文忽略了后效期火藥氣體對彈丸的加速作用，假設(shè)其以恒定速度運動，因此彈丸位置均有所滯后。

圖7 彈丸出膛前不同時刻流場溫度分布Fig.7 Temperature contours in the flow field at different times begore the bullet leaves the muzzle

圖8 計算陰影圖與實驗陰影照片比較Fig.8 Computational shadowgraph compared with experimental shadowgraph

圖9 t=415 μs時刻的密度、壓力云圖Fig.9 Density and pressure contours at t=415 μs

彈丸出膛后，負(fù)氧平衡的火藥氣體迅速射入初始流場，與環(huán)境中的氧氣接觸，在一定的溫度、壓力條件下可能發(fā)生化學(xué)反應(yīng)，甚至是劇烈燃燒。本文采用溫度T以及OH的質(zhì)量分?jǐn)?shù)w(OH)分布表征化學(xué)反應(yīng)劇烈程度，彈丸出膛后不同時刻分布云圖如圖10所示，可見算法成功捕捉到初始沖擊波、膛口沖擊波、入射斜激波、馬赫盤、彈底激波等復(fù)雜波系，清晰地展現(xiàn)了激波的傳播規(guī)律以及膛口焰的發(fā)展過程。圖11為彈丸與膛口間軸線溫度分布曲線。

圖10 不同時刻溫度(上)和OH質(zhì)量分?jǐn)?shù)(下)分布云圖Fig.10 Temperature contours (top) and mass fraction contours f OH (bottom) at different times

圖11 溫度沿軸線的分布曲線Fig.11 Temperature distribution along the axis

圖12 t=400 μs時刻N2質(zhì)量分?jǐn)?shù)分布及流線Fig.12 N2 mass fraction contours and streamlines at t=400 μs

隨著彈丸運動，膛內(nèi)高溫高壓火藥氣體進入初始流場而急劇膨脹，其速度達到2 km/s，迅速超越和包圍彈丸，形成相交激波和彈底激波，并且與空氣接觸面形成相對光滑的初始火焰陣面，如圖10(a)所示。由于初始流場的引導(dǎo)作用，加速了火藥氣體的軸向運動，使得膛口沖擊波呈冠狀，并迅速趕超初始沖擊波。由于邊界失穩(wěn)，氧氣被卷入，火焰陣面開始扭曲，如圖10(b)所示。圖10(c)中包含相交激波、馬赫盤在內(nèi)的完整欠膨脹射流結(jié)構(gòu)形成，并且膛口沖擊波基本湮沒初始流場。隨后彈丸穿越膛口沖擊波的過程中,如圖10(d)～10(e)所示，馬赫盤的生長不再受其影響，射流自由膨脹。t=320 μs時刻，由于自由膨脹，溫度大幅下降的火藥氣體經(jīng)過馬赫盤的再壓縮，溫度顯著上升，達到1 600 K,如圖11所示，但是由于馬赫盤下游區(qū)域缺乏氧氣，因此該區(qū)域并無劇烈化學(xué)反應(yīng)，但是由于高溫火藥氣體輻射出可見光，可形成中間焰。三波點附近出現(xiàn)向上游回旋的主渦環(huán)，將氧氣卷入射流邊界內(nèi)，并發(fā)生劇烈化學(xué)反應(yīng)，溫度達1 900 K，形成二次焰，此外可見OH質(zhì)量分?jǐn)?shù)分布并不均勻，而是呈卷曲的線狀，說明火藥氣體與環(huán)境的質(zhì)量交換仍不充分，化學(xué)反應(yīng)仍局限在有限的接觸面上。當(dāng)彈丸穿越膛口沖擊波后，對膛口流場影響進一步削弱，流場中高溫區(qū)域顯著增多，峰值達2 000 K以上，如圖10(f)～10(h)所示。此外受彈丸尾流的影響，部分馬赫盤后低速的火藥氣體被帶入下游，與環(huán)境中氧氣接觸，發(fā)生化學(xué)反應(yīng)。

N2在本文所采用的化學(xué)反應(yīng)機理中為惰性氣體，因此其質(zhì)量分?jǐn)?shù)分布可一定程度反應(yīng)火藥氣體與空氣質(zhì)量交換的強弱，t=400 μs時刻N2質(zhì)量分?jǐn)?shù)分布及流線如圖12所示，可見由于邊界不穩(wěn)定性形成了多個渦環(huán)，在其的作用下，射流邊界以及馬赫盤下游區(qū)域部分空氣被卷吸進入火藥氣體內(nèi)部，除三波點附近的主渦環(huán)外，馬赫盤下游、射流邊界上游存在著一系列小的渦環(huán)，這些渦環(huán)不斷生長、破碎，這些都使得原本清晰、簡單的接觸面卷曲，進而使火焰陣面變得更復(fù)雜。

6 結(jié) 論

改進了局部點云重構(gòu)方法，保證了點云均勻性，提高了計算效率，并建立了可用于包含大位移運動邊界，非平衡化學(xué)反應(yīng)流場數(shù)值模擬的無網(wǎng)格算法。通過對圓柱繞流以及激波誘導(dǎo)燃燒流場的計算驗證了算法的正確性。最后對12.7 mm高射機槍膛口二次燃燒流場進行了數(shù)值模擬，其結(jié)果清晰地展現(xiàn)了彈丸射出過程中，膛口初始流場、火藥氣體射流流場同高速運動彈丸的強烈耦合與相互作用；完整地再現(xiàn)了膛口流場復(fù)雜波系結(jié)構(gòu)變化的動力學(xué)過程以及二次焰分布的時空特征。

本文算法為包含任意位移動邊界化學(xué)反應(yīng)流的數(shù)值研究提供了有效的工具。

[1] Cai Xiao-wei, Tan Jun-jie, Ma Xin-jian, et al. Application of hybrid Cartesian grid and gridless approach to moving boundary flow problems[J]. International Journal for Numerical Method in Fluid, 2013,72(9):994-1013.

[2] 盛鳴劍,葉正寅,蔣超奇.附面層修正應(yīng)用于無網(wǎng)格算法的研究[J].計算力學(xué)學(xué)報,2011,28(6):920-925. Sheng Ming-jian, Ye Zheng-yin, Jiang Chao-qi. Study of boundary layer solution coupled with gridless method[J]. Chinese Journal of Computational Mechanics, 2011,28(6):920-925.

[3] Hashemi M Y, Jahangirian A. An efficient implicit mesh-less method for compressible flow calculations[J]. International Journal for Numerical Methods in Fluids, 2011,67(6):754-770.

[4] 周星.含動邊界復(fù)雜非定常流動的無網(wǎng)格算法研究[D].南京:南京理工大學(xué),2012.

[5] 馬新建.最下二乘無網(wǎng)格及其重疊點云法在CFD中的應(yīng)用研究[D].南京:南京理工大學(xué),2012.

[6] Kirshman D J, Liu F. A gridless boundary condition method for the solution of the Euler equations on embedded Cartesian meshes with multigrid[J]. Journal of Computational Physics, 2004,201(1):119-147.

[7] Luo H, Baum J D, L?hner R. A hybrid Cartesian grid and gridless method for compressible flows[J]. Journal of Computational Physics, 2006,214(2):618-632.

[8] 吳偉,許厚謙.激波誘導(dǎo)燃燒流場模擬的無網(wǎng)格算法[J].力學(xué)與實踐,2013,35(6):19-23. Wu Wei, Xu Hou-qian. Meshless method for numerical simulation of shock-induced combustion[J]. Mechanics in Engineering, 2013,35(6):19-23.

[9] Katz A, Jameson A. A comparison of various meshless schemes within a unified algorithm[R]. AIAA Paper, 2009-596,2009.

[10] 代淑蘭.復(fù)雜化學(xué)反應(yīng)流并行數(shù)值模擬[D].南京:南京理工大學(xué),2008.

[11] Soetrisno M, Imlay S T. Simulation of the flow field of a ram accelerator[R]. AIAA Paper, 1991-1915,1991.

(責(zé)任編輯 張凌云)

Numerical simulation of a muzzle flow field involving chemical reactions based on gridless method

Wu Wei, Xu Hou-qian, Wang Liang, Xue Rui

(SchoolofPowerEngineering,NanjingUniversityofScienceandTechnology,Nanjing210094,Jiangsu,China)

A gridless method for simulation of reactive flows involving moving boundaries was investigated based on the linear basis least-squares gridless method. The Euler equations of arbitrary Lagrangian-Eulerian form were employed as governing equations. The numerical flux and chemical sources were calculated by the multi-component HLLC ((Harten-Lax-van Leer-Contact) scheme and finite rate reaction model, respectively. An elevated restructuring technique of local cloud was adopted to deal with the moving boundaries. The front advance method of fictitious boundaries was used during the restructuring. The flow around the cylinder and the shock-induced combustion flow field were simulated to validate the accuracy firstly. The muzzle flow of a 12.7 mm machine gun was simulated. The computational shadowgraphs agree well with the experimental photographs, the density and pressure contours are in agreement with the results by the unstructured mesh method. The numerical results show the coupling and interaction progress in the initial muzzle flow field, the under-expanding jet of explosive gas and high-speed projectile clearly, and the temporal, spatial distribution characteristics of the muzzle flash are reappeared plainly.

mechanics of explosion; dynamic cloud; gridless method; muzzle flow field; chemical non-equilibrium flow

10.11883/1001-1455(2015)05-0625-08

2014-02-19；

2014-06-27

吳 偉(1990— )，男，博士研究生,wuwei.njust.804@gmail.com。

O381 國標(biāo)學(xué)科代碼： 13035

A