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

（南京理工大學(xué) 能源與動(dòng)力工程學(xué)院，南京210094）

傳統(tǒng)的數(shù)值方法主要包括有限元、有限差分、有限體積等，上述方法原理各不相同，但都需要采用網(wǎng)格劃分求解域，因而其應(yīng)用受到了網(wǎng)格的制約，如結(jié)構(gòu)網(wǎng)格處理相對(duì)復(fù)雜外形的求解區(qū)域能力有限；非結(jié)構(gòu)網(wǎng)格雖然對(duì)于復(fù)雜外形的處理靈活方便，但對(duì)于變形、動(dòng)邊界及不連續(xù)等復(fù)雜情形，生成高質(zhì)量的網(wǎng)格可能花費(fèi)與計(jì)算相當(dāng)?shù)臅r(shí)間；而笛卡爾網(wǎng)格雖然方便、快速，但是對(duì)于邊界的處理卻較為復(fù)雜、低效。為了擺脫網(wǎng)格的束縛，Batina［1］首先將無(wú)網(wǎng)格方法引入到計(jì)算流體力學(xué)領(lǐng)域，采用一系列點(diǎn)云離散求解域。無(wú)網(wǎng)格方法無(wú)需將節(jié)點(diǎn)連成網(wǎng)格，布點(diǎn)快捷，無(wú)需人為干涉，對(duì)于復(fù)雜外形更具靈活性。近十幾年，無(wú)網(wǎng)格方法受到國(guó)內(nèi)外學(xué)者的關(guān)注、研究，其內(nèi)容主要包含自動(dòng)布點(diǎn)生成點(diǎn)云，發(fā)展適用無(wú)網(wǎng)格方法的高激波分辨率、高精度數(shù)值格式，處理非定常問題的無(wú)網(wǎng)格方法以及混合網(wǎng)格方法等［2－5］，然而，現(xiàn)有的無(wú)網(wǎng)格方法多見應(yīng)用于空氣動(dòng)力學(xué)領(lǐng)域，對(duì)于航空航天、兵器等領(lǐng)域中廣泛存在的復(fù)雜化學(xué)反應(yīng)流場(chǎng)的模擬還顯得無(wú)能為力，對(duì)上述問題展開深入研究具有重要的意義和應(yīng)用價(jià)值。

本文基于PC機(jī)集群，研究了耦合化學(xué)反應(yīng)計(jì)算的并行無(wú)網(wǎng)格方法，用于包含復(fù)雜激波、燃燒波的非平衡化學(xué)反應(yīng)流的數(shù)值模擬。其基本變量解函數(shù)采用線性函數(shù)描述，空間導(dǎo)數(shù)采用最小二乘擬合逼近，將多組分HLLC格式引入無(wú)網(wǎng)格方法用于數(shù)值通量的計(jì)算，時(shí)間項(xiàng)采用4階Runge－Kutta法顯式推進(jìn)，反應(yīng)源項(xiàng)采用有限速率反應(yīng)模型處理，進(jìn)程間的數(shù)據(jù)通信采用MPI。采用上述無(wú)網(wǎng)格方法對(duì)激波誘導(dǎo)燃燒流場(chǎng)以及高速運(yùn)動(dòng)彈丸誘導(dǎo)斜爆轟流場(chǎng)進(jìn)行了數(shù)值模擬。

1 控制方程

本文采用含化學(xué)反應(yīng)源項(xiàng)的二維多組分Euler方程，具體形式如下：

守恒變量U，對(duì)流通量F、G，化學(xué)反應(yīng)源項(xiàng)W，軸對(duì)稱源項(xiàng)S的定義如下：

式中：t為流動(dòng)時(shí)間；ρs為組分s的質(zhì)量密度；N為組分總數(shù)；ρ為混合氣體的質(zhì)量密度；u，v分別為x，y方向速度分量；p為混合氣體總壓；ωs為化學(xué)反應(yīng)引起組分s的質(zhì)量生成率；ρE為單位體積總能。

式中：hs為組分s的焓值。各組分的焓值是溫度的函數(shù)，可按下多項(xiàng)式擬合獲得：

式中：Ru為通用氣體常數(shù)；Ms為組分s的摩爾質(zhì)量；T為混合氣體溫度；系數(shù)a1～a6可由JANAF表［6］查得。反應(yīng)中認(rèn)為混合氣體及各組分氣體滿足理想氣體狀態(tài)方程：

式中：M為混合氣體的平均摩爾質(zhì)量。

2 數(shù)值方法

2.1 對(duì)流通量計(jì)算

本文假設(shè)流動(dòng)解函數(shù)滿足線性關(guān)系，通過最小二乘擬合逼近空間導(dǎo)數(shù)。由此可得任意中心點(diǎn)i的對(duì)流通量為

形函數(shù)bij、cij以及空間導(dǎo)數(shù)具體求解過程參見文獻(xiàn)［5］，其中，中心點(diǎn)i及其衛(wèi)星點(diǎn)j的中點(diǎn)的數(shù)值通量Wij采用多組分HLLC格式計(jì)算，具體如下所示：

式中：下標(biāo)K根據(jù)式（7）取i或j，Si、Sj分別為中心點(diǎn)i、中心點(diǎn)j的波速，SM為接觸間斷面速度，lx、ly分別為單位法向量在x、y方向的分量，unK為中心點(diǎn)i或j處流體速度的法向分量，p＊按下式計(jì)算：

具體計(jì)算參見文獻(xiàn)［4－5］，在此不再贅述。

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

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

式中：ν′sm、ν″sm為組分s在反應(yīng)m中的化學(xué)反應(yīng)當(dāng)量系數(shù)，χs為組分s的化學(xué)符號(hào)，正向反應(yīng)速率常數(shù)Kfm可由Arrhenius公式計(jì)算獲得：

式中：Am為指前因子，bm為溫度因子，Em為活化能。逆向反應(yīng)速率常數(shù)Kbm可由反應(yīng)m的平衡常數(shù)或Arrhenius公式計(jì)算獲得。任意組分s的化學(xué)反應(yīng)源項(xiàng)可由下式計(jì)算：

式中：NR為化學(xué)反應(yīng)總數(shù)。

2.3 時(shí)間離散

對(duì)于任意中心點(diǎn)i，式（1）可寫成如下半離散形式：

式中：R為點(diǎn)i的殘差矢量。對(duì)式（14）可采用4階Runge－Kutta法進(jìn)行顯式時(shí)間推進(jìn)，形式如下：

式中：l＝1，2，3，4，φl(shuí)分別為1／4、1／3、1／2、1。為了保證計(jì)算的穩(wěn)定，采用各進(jìn)程中最小時(shí)間步長(zhǎng)作為全局時(shí)間步長(zhǎng)Δt。流動(dòng)與化學(xué)反應(yīng)解耦計(jì)算，為克服“剛性”問題，將化學(xué)反應(yīng)時(shí)間步長(zhǎng)作進(jìn)一步細(xì)分處理。

3 并行處理

本文并行計(jì)算在多臺(tái)通過D－Link集線器互連的Dell臺(tái)式機(jī)上進(jìn)行，進(jìn)程之間的數(shù)據(jù)通信均采用具有良好可移植性、可擴(kuò)展性的MPI。并行計(jì)算中負(fù)載平衡是影響計(jì)算效率的重要因素。合理的分配策略是保證并行效率的關(guān)鍵。為了保證各進(jìn)程分配均勻，本文采用相對(duì)簡(jiǎn)單的一維劃分，將離散點(diǎn)進(jìn)行分區(qū)編號(hào)，根據(jù)進(jìn)程數(shù)平均分配，各進(jìn)程保存各自分區(qū)的離散點(diǎn)，并生成相應(yīng)的分區(qū)通信邊界，建立向相鄰進(jìn)程發(fā)送、接受數(shù)據(jù)離散點(diǎn)鏈表。計(jì)算每推進(jìn)1個(gè)時(shí)間步長(zhǎng)，各進(jìn)程向相鄰進(jìn)程發(fā)送并從相鄰進(jìn)程接受通信邊界節(jié)點(diǎn)數(shù)據(jù)。

4 算例

共選取2個(gè)算例驗(yàn)證本文算法的有效性，分別為激波誘導(dǎo)燃燒和高速飛行彈丸誘導(dǎo)爆轟實(shí)驗(yàn)，算例1是半徑R為7.5mm的球－柱體鈍頭彈射入等化學(xué)當(dāng)量比的H2／Air預(yù)混氣體，算例2為尖頭彈（簡(jiǎn)化結(jié)構(gòu)如圖1所示）射入等化學(xué)當(dāng)量比的CH4／Air預(yù)混氣體，壓力p0、溫度T0、速度v0如表1所示。氫氣／空氣采用7組分8反應(yīng)模型，正／逆反應(yīng)速率常數(shù)計(jì)算系數(shù)取自文獻(xiàn)［7］，甲烷／空氣采用13組分19反應(yīng)模型，正向反應(yīng)速率常數(shù)計(jì)算系數(shù)同文獻(xiàn)［8］，逆向反應(yīng)速率由平衡常數(shù)計(jì)算。本文計(jì)算離散點(diǎn)通過非結(jié)構(gòu)網(wǎng)格無(wú)網(wǎng)格化方法獲得。

圖1 駐定斜爆轟波實(shí)驗(yàn)彈丸結(jié)構(gòu)

表1 實(shí)驗(yàn)條件

算例1采用非均勻分布的121 226個(gè)節(jié)點(diǎn)離散流場(chǎng)，其中駐點(diǎn)流線上布點(diǎn)140個(gè)。圖2、圖3分別為采用4進(jìn)程計(jì)算得到的等溫線和H2O質(zhì)量分?jǐn)?shù)分布，圖2中“■”為實(shí)驗(yàn)照片中激波陣面的位置，可見自由來(lái)流速度大于臨界爆轟速度，激波陣面與燃燒陣面耦合。圖4、圖5為沿駐點(diǎn)流線無(wú)量綱壓力和溫度的分布曲線，可見激波后溫度、壓力迅速上升并達(dá)到平衡狀態(tài)，同文獻(xiàn)［9］采用網(wǎng)格方法得到的結(jié)果吻合較好。此外本文算法成功捕捉到壓力的VonNeumann尖峰。

圖2 等溫線分布

圖3 H2O質(zhì)量分?jǐn)?shù)分布

圖4 沿駐點(diǎn)流線無(wú)量綱壓力分布

圖5 沿駐點(diǎn)流線無(wú)量綱溫度分布

算例2采用非均勻分布的207 532個(gè)節(jié)點(diǎn)離散，對(duì)彈丸表面附近的節(jié)點(diǎn)進(jìn)行了加密處理。圖6為采用3進(jìn)程計(jì)算得到的密度云圖，作為比較，同樣給出了實(shí)驗(yàn)陰影照片。

圖7上下分別為無(wú)網(wǎng)格方法與非結(jié)構(gòu)網(wǎng)格［8］方法計(jì)算得到的溫度云圖，可見彈丸頭部激波較弱，不足以點(diǎn)燃可燃混合氣體，僅形成駐定的斜激波，計(jì)算得到激波角為35.4°，同實(shí)驗(yàn)獲得的37.8°以及非結(jié)構(gòu)網(wǎng)格計(jì)算得到的36.5°吻合較好。圖8為CH4以及CO2質(zhì)量分?jǐn)?shù)分布云圖，可見僅在彈丸中部凹槽內(nèi)以及彈尾發(fā)生了化學(xué)反應(yīng)。

圖6 密度云圖與實(shí)驗(yàn)陰影照片對(duì)比

圖7 溫度云圖對(duì)比

圖8 CH4和CO2質(zhì)量分?jǐn)?shù)分布

5 效率分析

本文分別使用2～8個(gè)進(jìn)程對(duì)上述2個(gè)算例進(jìn)行了并行計(jì)算。由于并行程序同串行程序存在差異，將兩者進(jìn)行比較不夠合理，因而本文采用2進(jìn)程計(jì)算耗時(shí)為參考量計(jì)算加速比δ、并行效率η，其部分結(jié)果見表2。可見隨著進(jìn)程個(gè)數(shù)NP的增加，并行效率逐漸下降，這是由于NP的增加導(dǎo)致各進(jìn)程分配有效計(jì)算節(jié)點(diǎn)的減少，使得計(jì)算與通信的耗時(shí)比減小。此外雖然算例2離散點(diǎn)個(gè)數(shù)多于算例1，但是需要計(jì)算化學(xué)反應(yīng)的節(jié)點(diǎn)卻遠(yuǎn)少于算例1，其計(jì)算時(shí)間與通信時(shí)間相比所占比重較低，因此算例2的并行效率均低于算例1。

表2 并行加速比及并行效率

6 結(jié)束語(yǔ)

本文基于MPI并行壞境，發(fā)展了非平衡化學(xué)反應(yīng)流數(shù)值模擬的并行無(wú)網(wǎng)格算法，并對(duì)激波誘導(dǎo)燃燒和高速飛行彈丸誘導(dǎo)爆轟流場(chǎng)進(jìn)行了數(shù)值模擬，其結(jié)果同實(shí)驗(yàn)以及網(wǎng)格方法結(jié)果吻合較好，并且獲得了較為理想的并行效率，表明了計(jì)算形式簡(jiǎn)單的多組分HLLC格式具有較高的激波分辨率，驗(yàn)證了本文并行無(wú)網(wǎng)格算法的正確性。本文研究不僅拓展了無(wú)網(wǎng)格方法處理復(fù)雜流場(chǎng)的能力，豐富了非平衡反應(yīng)流模擬的數(shù)值方法，而且為航空航天、兵器等領(lǐng)域中眾多化學(xué)反應(yīng)流的大規(guī)模數(shù)值計(jì)算提供了有效的工具。

［1］BATINA J T.A gridless Euler／Navier－Stokes solution algorithm for complex－aircraft applications，AIAA 93－0333［R］.1993.

［2］孫迎丹，王剛，葉正寅.AUSM＋－up格式在無(wú)網(wǎng)格算法中的推廣［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2005，23（4）：511－515.SUN Ying－dan，WANG Gang，YE Zheng－yin.Extending the AUSM＋－up scheme into the gridless method［J］.Acta Aerodynamica Sinica，2005，23（4）：511－515.（in Chinese）

［3］WANG H，CHEN H Q，PERIAUX J.A study of gridless method with dynamic clouds of points for solving unsteady CFD problems in aerodynamics［J］.Int J Meth Fluids，2010，64：98－118.

［4］周星，許厚謙.無(wú)網(wǎng)格算法在膛口流場(chǎng)數(shù)值模擬中的應(yīng)用［J］.彈道學(xué)報(bào)，2011，23（1）：23－26.ZHOU Xing，XU Hou－qian.Application of gridless method in numerical simulation of muzzle flow field［J］.Journal of ballistics，2011，23（1）：23－26.（in Chinese）

［5］ZHOU X，XU H Q.Gridless method for unsteady flows involving moving discrete points and its applications［J］.Engineering Applications of Computational Fluid Mechanics，2010，4（2）：276－286.

［6］STULL D R，PROPHET H.JANAF thermochemical tables［J］.Journal of Physical and Chemical Reference Data，1975，4（1）：1－176.

［7］EVANS J S，SCHEXNAYDER C J.Influence of chemical kinetics and unmixedness on burning in supersonic hydrogen flames［J］.AIAA，1980，18（2）：188－193.

［8］代淑蘭，許厚謙.駐定斜爆轟波并行數(shù)值模擬［J］.爆炸與沖擊，2006，26（6）：572－576.DAI Shu－lan，XU Hou－qian.Numerical simulation of standing oblique detonation wave using parallel computation method［J］.Explosion and Shock Waves，2006，26（6）：572－576.（in Chinese）

［9］孫明波，梁劍寒，王振國(guó).非平衡流解耦方法及其計(jì)算激波誘導(dǎo)燃燒的應(yīng)用驗(yàn)證［J］.航空動(dòng)力學(xué)報(bào)，2008，23（11）：2 055－2 061.SUN Ming－bo，LIANG Jian－h(huán)an，WANG Zhen－guo.Validation of an uncoupled solver of non－equilibrium flow for shock－induced combustion［J］.Journal of Aerospace Power，2008，23（11）：2 055－2 061.（in Chinese）