化學反應流模擬的并行無網格方法

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

（南京理工大學 能源與動力工程學院，南京210094）

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

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

1 控制方程

本文采用含化學反應源項的二維多組分Euler方程，具體形式如下：

守恒變量U，對流通量F、G，化學反應源項W，軸對稱源項S的定義如下：

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

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

式中：Ru為通用氣體常數；Ms為組分s的摩爾質量；T為混合氣體溫度；系數a1～a6可由JANAF表［6］查得。反應中認為混合氣體及各組分氣體滿足理想氣體狀態方程：

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

2 數值方法

2.1 對流通量計算

本文假設流動解函數滿足線性關系，通過最小二乘擬合逼近空間導數。由此可得任意中心點i的對流通量為

形函數bij、cij以及空間導數具體求解過程參見文獻［5］，其中，中心點i及其衛星點j的中點的數值通量Wij采用多組分HLLC格式計算，具體如下所示：

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

具體計算參見文獻［4－5］，在此不再贅述。

2.2 化學反應模型

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

式中：ν′sm、ν″sm為組分s在反應m中的化學反應當量系數，χs為組分s的化學符號，正向反應速率常數Kfm可由Arrhenius公式計算獲得：

式中：Am為指前因子，bm為溫度因子，Em為活化能。逆向反應速率常數Kbm可由反應m的平衡常數或Arrhenius公式計算獲得。任意組分s的化學反應源項可由下式計算：

式中：NR為化學反應總數。

2.3 時間離散

對于任意中心點i，式（1）可寫成如下半離散形式：

式中：R為點i的殘差矢量。對式（14）可采用4階Runge－Kutta法進行顯式時間推進，形式如下：

式中：l＝1，2，3，4，φl分別為1／4、1／3、1／2、1。為了保證計算的穩定，采用各進程中最小時間步長作為全局時間步長Δt。流動與化學反應解耦計算，為克服“剛性”問題，將化學反應時間步長作進一步細分處理。

3 并行處理

本文并行計算在多臺通過D－Link集線器互連的Dell臺式機上進行，進程之間的數據通信均采用具有良好可移植性、可擴展性的MPI。并行計算中負載平衡是影響計算效率的重要因素。合理的分配策略是保證并行效率的關鍵。為了保證各進程分配均勻，本文采用相對簡單的一維劃分，將離散點進行分區編號，根據進程數平均分配，各進程保存各自分區的離散點，并生成相應的分區通信邊界，建立向相鄰進程發送、接受數據離散點鏈表。計算每推進1個時間步長，各進程向相鄰進程發送并從相鄰進程接受通信邊界節點數據。

4 算例

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

圖1 駐定斜爆轟波實驗彈丸結構

表1 實驗條件

算例1采用非均勻分布的121 226個節點離散流場，其中駐點流線上布點140個。圖2、圖3分別為采用4進程計算得到的等溫線和H2O質量分數分布，圖2中“■”為實驗照片中激波陣面的位置，可見自由來流速度大于臨界爆轟速度，激波陣面與燃燒陣面耦合。圖4、圖5為沿駐點流線無量綱壓力和溫度的分布曲線，可見激波后溫度、壓力迅速上升并達到平衡狀態，同文獻［9］采用網格方法得到的結果吻合較好。此外本文算法成功捕捉到壓力的VonNeumann尖峰。

圖2 等溫線分布

圖3 H2O質量分數分布

圖4 沿駐點流線無量綱壓力分布

圖5 沿駐點流線無量綱溫度分布

算例2采用非均勻分布的207 532個節點離散，對彈丸表面附近的節點進行了加密處理。圖6為采用3進程計算得到的密度云圖，作為比較，同樣給出了實驗陰影照片。

圖7上下分別為無網格方法與非結構網格［8］方法計算得到的溫度云圖，可見彈丸頭部激波較弱，不足以點燃可燃混合氣體，僅形成駐定的斜激波，計算得到激波角為35.4°，同實驗獲得的37.8°以及非結構網格計算得到的36.5°吻合較好。圖8為CH4以及CO2質量分數分布云圖，可見僅在彈丸中部凹槽內以及彈尾發生了化學反應。

圖6 密度云圖與實驗陰影照片對比

圖7 溫度云圖對比

圖8 CH4和CO2質量分數分布

5 效率分析

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

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

6 結束語

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

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

［2］孫迎丹，王剛，葉正寅.AUSM＋－up格式在無網格算法中的推廣［J］.空氣動力學學報，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］周星，許厚謙.無網格算法在膛口流場數值模擬中的應用［J］.彈道學報，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］代淑蘭，許厚謙.駐定斜爆轟波并行數值模擬［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］孫明波，梁劍寒，王振國.非平衡流解耦方法及其計算激波誘導燃燒的應用驗證［J］.航空動力學報，2008，23（11）：2 055－2 061.SUN Ming－bo，LIANG Jian－han，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）