含化學反應膛口流場的無網格數值模擬*

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

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

?

含化學反應膛口流場的無網格數值模擬*

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

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

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

爆炸力學；動態點云；無網格方法；膛口流場；非平衡反應流

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

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

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

1 控制方程

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

(1)

式中：U為守恒變量，F、G為對流通量，W為化學反應源項，S為軸對稱源項。

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

2 點云重構

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

3 數值方法

3.1 空間離散

本文中發展的無網格方法采用線性基函數最小二乘擬合[4,9]直接求解各離散點的空間導數。假設流動基本變量滿足如下線性關系：

(2)

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

(3)

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

(4)

進而可得中心點i的空間導數為：

(5)

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

(6)

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

(7)

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

(9)

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

(10)

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

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

(11)

(12)

(13)

(14)

3.2 化學反應模型

化學反應源項采用有限速率反應模型，反應體系中的任意反應可表示為：

(15)

(16)

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

(17)

式中：NR為化學反應總數，Mi為組分i的摩爾質量。

3.3 時間項及邊界條件

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

4 驗證算例

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

4.1 圓柱繞流流場

為了驗證點云重構技術的有效性，對圓柱繞流問題進行了模擬，計算域為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，計算中凍結化學反應。t=0.8 ms時的壓力云圖如圖1所示，圓柱上表面的量綱一壓力p/p0分布如圖2所示，可見采用動態點云重構計算得到的結果同圓柱靜止時結果吻合較好。

圖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 激波誘導燃燒流場

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

圖3 溫度云圖Fig.3 Temperature contours

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

圖5 駐點流線上組分質量分數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的高射機槍膛口流場進行了數值模擬，膛管的內徑為13.7 mm，外徑為31.0 mm，膛管長為1.08 m，計算中對彈丸結構進行了適當的簡化，其出膛速度為810 m/s。外流場取0.32 m×0.8 m矩形區域，流場共布點228 247個，某時刻流場局部區域離散點分布見圖6。定義彈底到達膛口瞬間為0時刻。僅在彈丸出膛后開啟化學反應計算，其化學反應采用H2-CO-O29組分11步反應機理[10]。

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

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

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

vx=v0x/L

(18)

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

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

圖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

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

圖10 不同時刻溫度(上)和OH質量分數(下)分布云圖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質量分數分布及流線Fig.12 N2 mass fraction contours and streamlines at t=400 μs

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

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

6 結 論

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

本文算法為包含任意位移動邊界化學反應流的數值研究提供了有效的工具。

[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] 盛鳴劍,葉正寅,蔣超奇.附面層修正應用于無網格算法的研究[J].計算力學學報,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] 周星.含動邊界復雜非定常流動的無網格算法研究[D].南京:南京理工大學,2012.

[5] 馬新建.最下二乘無網格及其重疊點云法在CFD中的應用研究[D].南京:南京理工大學,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] 吳偉,許厚謙.激波誘導燃燒流場模擬的無網格算法[J].力學與實踐,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] 代淑蘭.復雜化學反應流并行數值模擬[D].南京:南京理工大學,2008.

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

(責任編輯 張凌云)

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 國標學科代碼： 13035

A