楊浩森，張 斌，劉 洪，陳 方

（上海交通大學(xué) 航空航天學(xué)院，上海 200240）

0 引 言

直接模擬蒙特卡洛法（DSMC）的基本原理是采用有限個模擬分子代替真實氣體分子，直接模擬微觀態(tài)的氣體分子運動和碰撞，是一種更接近于氣體真實的運動狀態(tài)和能量交換的模擬方法。鑒于這種直接的物理模型，與基于連續(xù)介質(zhì)假設(shè)的CFD方法相比，DSMC在模擬稀薄流區(qū)流動和涉及熱化學(xué)非平衡的流動時具有天然的優(yōu)勢。

由于DSMC以統(tǒng)計學(xué)為基礎(chǔ)，則作為流場統(tǒng)計母體的模擬分子數(shù)量對流場最后的統(tǒng)計結(jié)果將產(chǎn)生直接的影響。這種直接的影響主要表現(xiàn)在計算效率和計算精度兩個方面。概括來講，就是在保證計算精度的同時盡量減少流場中使用的模擬分子數(shù)以降低計算開銷。

對于不涉及熱化學(xué)非平衡現(xiàn)象的流場，流場中模擬分子的主要行為是相互碰撞以實現(xiàn)速度的重新分布，這種行為只涉及分子平動動能，物理模型相對簡單。Bird和Alexander等人對此類流動中模擬分子數(shù)對DSMC計算的影響做了定性分析，建議每個流場網(wǎng)格單元中布置15－25個模擬分子以得到合理的流場計算結(jié)果［1－3］，該研究結(jié)論作為此類流場模擬的推薦值被廣泛使用，Pei－Yuan Tzeng等人在此基礎(chǔ)上以Rayleigh－Bénard流動為例分析了模擬分子數(shù)對流場渦結(jié)構(gòu)準(zhǔn)確性的影響［4］。2009年，Bird對全新的DSMC07程序中網(wǎng)格模擬分子數(shù)做了更為準(zhǔn)確全面地分析，得出從1維到3維流場理論計算“改進因子”與網(wǎng)格模擬分子數(shù)的關(guān)系，結(jié)果表明在2維不含熱化學(xué)非平衡流場的DSMC計算中，相同網(wǎng)格模擬分子數(shù)對計算的改進效果并不會因為流場物理條件的不同而不同。結(jié)合計算效率，Bird建議，在二維不含熱化學(xué)非平衡的流場模擬中網(wǎng)格模擬分子數(shù)應(yīng)該在20－30個，此結(jié)論與流場物理條件無關(guān)［5］。

上述討論僅僅局限在不涉及熱化學(xué)非平衡的流場模擬中。但是，在研究例如高超聲速飛行器和飛行器再入等問題時，流場中會出現(xiàn)顯著的熱化學(xué)非平衡的現(xiàn)象，這類現(xiàn)象對流場特征有十分重要的影響［6－7］。與傳統(tǒng)方法不同，DSMC在模擬此類非平衡現(xiàn)象時，要求模擬分子除了實現(xiàn)彈性碰撞之外，還要準(zhǔn)確模擬分子內(nèi)能的激發(fā)與松弛過程。概括而言，與普通流場中DSMC主要考慮平動能不同，模擬分子的碰撞行為將同時涉及平動，轉(zhuǎn)動和振動三種能量的相互傳遞。對于此類更加復(fù)雜的流場模擬，物理條件對分子內(nèi)能的影響是顯著的，要求在不同的物理狀態(tài)下（主要是溫度），碰撞模擬應(yīng)該符合不同的能量松弛規(guī)律，不能如普通流場一概視之［8］。

本文使用DSMC模擬了空氣中最易發(fā)生化學(xué)反應(yīng)的三種氮氧分子：氮氣，氧氣和一氧化氮在熱浴條件下的離解反應(yīng)，并以此來研究在涉及熱化學(xué)非平衡的流場中模擬分子數(shù)對計算結(jié)果的影響。并從分子能態(tài)分布和能量松弛等方面分析了這種影響產(chǎn)生的內(nèi)在原因。

1 DSMC熱化學(xué)非平衡模擬方法

1.1 概述

根據(jù)分子動力學(xué)和化學(xué)動力學(xué)的基本原理，DSMC模擬熱化學(xué)非平衡主要分為兩個部分：熱力學(xué)能量傳遞和各類化學(xué)反應(yīng)［9］。這兩個過程均通過模擬分子之間的碰撞實現(xiàn)。在不涉及熱化學(xué)非平衡流場的計算中，模擬分子主要攜帶分子平動能，而與普通流場計算中的模擬分子相比，模擬熱化學(xué)非平衡流場要求每一個模擬分子具有相應(yīng)的內(nèi)能能態(tài)，具體如下式所示：

普通流場模擬分子與非平衡流場模擬分子的主要區(qū)別如表1所示。與普通流場不同，流場的非平衡現(xiàn)象主要就來源于分子的內(nèi)能能態(tài)的激發(fā)與松弛過程，對于平動能的傳能而言，一般情況下可以用彈性碰撞表示，而將涉及分子內(nèi)能傳遞的碰撞歸入非彈性碰撞類。

表1 模擬分子特性比較Table 1 Different between normal and nonequilibrium simulation molecules

上述所有的能量形式都存在相互獨立的能量松弛時間，在不考慮分子電離的情況下，三種主要的能量形式：平動能（Transitional Energy），轉(zhuǎn)動能（Rotational Energy）和振動能（Vibrational Energy）所對應(yīng)的松弛時間存在如下關(guān)系：

松弛時間的倒數(shù)被定義為松弛碰撞數(shù)，用以描述單位時間內(nèi)單個分子相應(yīng)能量松弛所需要的碰撞次數(shù)。一般而言，平動能只需要單次碰撞便可以實現(xiàn)松弛過程，即分子速度重新分布，此類碰撞模型相對簡單，而且平動能松弛與實際的物理條件無關(guān)。而其他兩種分子內(nèi)能能態(tài)卻相對比較復(fù)雜，在不同的流場溫度下，內(nèi)能松弛需要的時間是不同的［6］，亦即松弛碰撞數(shù)是不同的。因此，DSMC在不同的流場狀態(tài)下準(zhǔn)確模擬內(nèi)能松弛碰撞是影響流場終態(tài)準(zhǔn)確性的關(guān)鍵因素。

在進行分子各能態(tài)能量之間的相互有效傳遞之后，再根據(jù)碰撞總能和各類化學(xué)反應(yīng)的相應(yīng)特征判定化學(xué)反應(yīng)的進行。

1.2 DSMC模擬離解反應(yīng)

到目前為止，被廣泛使用的標(biāo)準(zhǔn)DSMC化學(xué)反應(yīng)能量模型是Bird所提出的TCE（Total Collision Energy）模型［10］。在這個模型中，離解氣體分子的離解反應(yīng)是由分子振動能主導(dǎo)的，根據(jù)Bird的TCE模型，對一個特定的分子碰撞對而言，碰撞能量分為分子平動動能，轉(zhuǎn)動能和振動能。在分子對碰撞的過程中，按照Larsen－Borgnakke模型進行各自由度之間的能量傳遞和再分配，最后的平衡態(tài)能量分布根據(jù)Boltzmann平衡分布給出［11－12］。

在DSMC模擬化學(xué)反應(yīng)的過程中，模擬分子之間的非彈性碰撞會導(dǎo)致分子各自由度能量的重新分配，根據(jù)化學(xué)反應(yīng)的細致平衡原理［9］，分子內(nèi)能（轉(zhuǎn)動能和振動能）在分子相對平動能的直接激發(fā)下對化學(xué)反應(yīng)做出相應(yīng)的貢獻。而對于離解反應(yīng)而言，由于上述能量之間的轉(zhuǎn)化，導(dǎo)致分子振動量子態(tài)從低能態(tài)向高能態(tài)躍遷，當(dāng)分子振動能態(tài)超出特定的離解能級時，才發(fā)生離解反應(yīng)［13］。本文分子振動能級分布采用式（3）所示的非簡諧振子模型，其中，相關(guān)的能級參數(shù)如表2所示。

表2 振動能級參數(shù)［14］Table 2 Parameters in vibrational energy model

對于離解反應(yīng)的模擬，一個重要的考量標(biāo)準(zhǔn)就是不同溫度下的平衡離解度α（Degree of Dissociation－DOD），它描述了在不同的平衡溫度下，某一類分子離解的程度，即離解原子質(zhì)量分?jǐn)?shù)。根據(jù)Bird和Hansen（1976）提供的數(shù)據(jù)［12］，這里不加推導(dǎo)的給出氮氣、氧氣和一氧化氮在不同溫度下的平衡離解度公式為：

2 計算結(jié)果及分析

2.1 DSMC程序驗證

首先對DSMC程序模擬三種氮氧分子離解反應(yīng)的正確性進行驗證。以標(biāo)準(zhǔn)狀態(tài)（1個大氣壓，0攝氏度）分子數(shù)密度為準(zhǔn)，為提高計算效率，采用一個封閉系統(tǒng)的熱浴反應(yīng)來模擬三種氣體在不同平衡溫度下的離解反應(yīng)。根據(jù)統(tǒng)計學(xué)原理，隨機系統(tǒng)的統(tǒng)計漲落僅依賴于總的模擬分子數(shù)，而與模擬分子的布置無關(guān)［15］，故將整個計算域視為一個網(wǎng)格，網(wǎng)格模擬分子數(shù)為7000個。全文DSMC程序參數(shù)如表3所示。

表3 本文全局算例計算參數(shù)Table 3 Global calculation parameters of DSMC

圖1所示的結(jié)果為本文使用的DSMC計算程序在表3中列舉的不同分子數(shù)密度下對三種氮氧分子進行熱浴的模擬結(jié)果。氮氣的模擬平衡溫度范圍大致為5000K～16000K，氧氣大致為2600K～9000K，一氧化氮大致為3000K～11000K，圖例中n0＝2.686×1023。由于分子的離解度與計算域的溫度直接相關(guān)，所以本文選擇簡單絕熱熱浴問題進行計算模擬，以盡可能多的排除其他影響因素從而突出模擬分子數(shù)的影響。

圖1 DSMC模擬氧氣和氮氣離解度曲線Fig.1 DSMC result of degree of dissociation

從圖1可以看出，由于三種分子化學(xué)鍵鍵能不同，則離解特征溫度存在差異。其中氮氣最為穩(wěn)定，整個離解溫度范圍最高；氧氣最易離解，而一氧化氮正好介于氮氣和氧氣之間。DSMC模擬結(jié)果在不同的分子數(shù)密度下和理論值符合得比較好，表明程序是可靠的。

進一步分析之前，分別對每一種氣體選擇一個標(biāo)準(zhǔn)算例作為判斷模擬分子數(shù)對計算結(jié)果影響的標(biāo)定。計算條件為：分子數(shù)密度為標(biāo)準(zhǔn)狀態(tài)下的分子數(shù)密度，模擬分子數(shù)為4000。圖2為氧氣，氮氣和一氧化氮標(biāo)準(zhǔn)算例的計算結(jié)果和理論值的誤差基本在10%以內(nèi)（圖中誤差限為5%）。

圖2 標(biāo)準(zhǔn)算例結(jié)果Fig.2 DSMC results of the standard examples

在圖2所示的標(biāo)準(zhǔn)算例中，與高溫高離解度的情況相比，低溫低離解度下DSMC模擬結(jié)果與理論值之間的誤差較大（計算結(jié)果顯示最大誤差達到了12%左右），其原因根據(jù)Gallis和Harvey的研究結(jié)果［16］，空氣中氣體的化學(xué)反應(yīng)并非由振動能唯一控制，在分子振動能還沒有被高度激發(fā)的狀態(tài)下，平動能和轉(zhuǎn)動能對化學(xué)反應(yīng)的影響將變得不可忽略，本文使用的DSMC程序模擬化學(xué)反應(yīng)是基于分子振動能的控制，故在低溫下模擬誤差較大。

2.2 模擬分子數(shù)對DSMC模擬結(jié)果的影響

圖3給出了在和標(biāo)準(zhǔn)算例相同的計算條件下，不同的模擬分子數(shù)在不同的平衡離解度下模擬的誤差（圖中的虛線表示上一節(jié)規(guī)定的相應(yīng)離解度下標(biāo)準(zhǔn)算例計算結(jié)果與理論值之間的誤差，以此作為計算結(jié)果被接受的標(biāo)定值）。

圖3中標(biāo)識出了在不同的平衡離解度下需要的模擬分子數(shù)下限，超過這個閾值之后，計算精度并沒有得到明顯的提升，但是由于模擬分子數(shù)增多帶來的計算開銷卻會相應(yīng)增大。綜合該離解度下的流場平衡溫度，總結(jié)得到的模擬粒子數(shù)需求如表4所示。

考慮到DSMC的統(tǒng)計誤差，故表3給出的分子數(shù)閾值以一個小范圍給出。很明顯，隨著溫度的升高，氮氧分子離解度上升，此時要求的模擬分子數(shù)也隨之升高。離解度為95%時要求的初始模擬分子數(shù)是5%離解度時的10倍左右，這一結(jié)果與Bird早期提出的對普通流場DSMC模擬分子數(shù)選取的統(tǒng)一建議標(biāo)準(zhǔn)不同。本文結(jié)果表明，化學(xué)反應(yīng)系統(tǒng)要求的模擬分子數(shù)根據(jù)物理條件不同應(yīng)當(dāng)作區(qū)別對待。

圖3 不同模擬分子數(shù)下離解度誤差Fig.3 Degree of disassociation error with different simulation molecules

表4 不同理論離解度下需要的模擬分子數(shù)Table 4 The effective number of simulation molecules in different DOD

（1）不同模擬分子數(shù)對分子能級分布的影響

根據(jù)Bird的理論，在計算穩(wěn)定之后，流場中的分子振動能級分布將滿足Boltzmann平衡分布，在理論離解度為75%的情況下三種氮氧氣體使用不同的模擬分子數(shù)最后流場分子達到的狀態(tài)如圖4所示。

圖4 三種分子在不同模擬分子數(shù)下的計算平衡能級分布Fig.4 Equilibrium energy distribution in different simulation molecules

根據(jù)上一節(jié)離解誤差的結(jié)論，理論離解度為75%時要求模擬分子數(shù)為900個。圖4顯示的結(jié)果表明，在使用900個模擬分子的情況下，最后DSMC分子能級抽樣統(tǒng)計結(jié)果和Boltzmann平衡分布一致。而在模擬分子數(shù)為150時，最后的分子能級分布偏離Boltzmann分布較遠，這一結(jié)果主要由兩個因素所致：其一，從統(tǒng)計學(xué)的角度出發(fā)，母體樣本數(shù)太少，導(dǎo)致的抽樣不準(zhǔn)確；其二，初始模擬分子數(shù)不足，導(dǎo)致分子碰撞傳能不夠充分，導(dǎo)致流場能量分布發(fā)生偏倚，這一點將在下一節(jié)做詳細討論。

（2）不同模擬分子數(shù)對有效碰撞率和振動松弛碰撞率的影響

分子之間的能量傳遞是通過分子對之間的碰撞實現(xiàn)的，根據(jù)Parker的理論［17］，在不同溫度下，分子振動能的松弛需要滿足不同的松弛碰撞數(shù)（即分子的振動能要完全松弛要求發(fā)生的平均碰撞次數(shù)），且在特定的溫度下，振動松弛碰撞數(shù)大致為一個確定值，溫度越高需要的松弛碰撞數(shù)越少。根據(jù)這個理論，本文定義了以下兩個考察參數(shù)，分別為：“有效碰撞率”（Effective Collision Rate）和“振動松弛碰撞率”（Vibrational Relaxation Collision Rate）。從模擬離解反應(yīng)的角度出發(fā)，在DSMC模擬過程中，將有雙原子分子參與的碰撞定義為“有效碰撞（NCmoleculeinvolved）”，而“有效碰撞率（νeffective）”描述的是“有效碰撞”占總碰撞數(shù)（NCtotal）的比例?！罢駝铀沙谂鲎猜剩é蛌ib）”則描述在“有效碰撞”中，發(fā)生能量松弛過程的碰撞數(shù)（NCvibrelax）比例，具體表述如下：

圖5給出了在平衡離解度分別為30%和75%兩種情況下三中氮氧分子碰撞傳能情況。從圖中可以明顯看出，分子有效碰撞率和分子振動松弛碰撞率隨著模擬分子數(shù)的增加，最后都各自趨于一個穩(wěn)定值，理論離解度為30%時大致穩(wěn)定的模擬分子數(shù)為300個，離解度75%時大致穩(wěn)定的模擬分子數(shù)為900個，這個結(jié)果和上文的結(jié)論一致。

模擬分子數(shù)過少，離解不足，流場中剩余的雙原子分子數(shù)量較多，此時有效碰撞率較高，另一方面，松弛分子碰撞對抽樣不足，不能達到足夠的松弛碰撞數(shù)，即振動松弛碰撞率較低。隨著模擬分子數(shù)的增加，計算離解度升高，此時剩余雙原子分子數(shù)量減少，則有效碰撞率下降。隨著模擬分子數(shù)的增加松弛碰撞抽樣增加，松弛碰撞率增大直至穩(wěn)定。這種趨勢表明，同一平衡狀態(tài)下，正確模擬的平衡碰撞率跟模擬分子數(shù)量無關(guān)，為一恒定值。這個結(jié)果與Parker的理論一致。

考慮極端的情況，在氣體幾乎完全離解的情況下，由于單原子并不具有內(nèi)部自由度，亦即不具有轉(zhuǎn)動和振動等內(nèi)能能態(tài)，從統(tǒng)計學(xué)的角度，如果離解度過高且沒有足夠的分子（雙原子分子）為程序提供建立平衡狀態(tài)所需要的采樣母體的規(guī)模，就會直接導(dǎo)致流場能量統(tǒng)計失準(zhǔn)，最直接的結(jié)果就是流場溫度的統(tǒng)計值不準(zhǔn)確。所以就要求起始時刻有足夠多的模擬分子基數(shù)，這樣，隨著離解反應(yīng)的進行，也能夠有足夠的未離解分子供程序進行采樣統(tǒng)計從而得出正確的流場狀態(tài)，而隨著模擬的進行，流場逐漸趨近穩(wěn)定狀態(tài)，這時即便是分子完全離解也不會影響最后的模擬結(jié)果。

圖5 不同模擬分子數(shù)下分子碰撞率Fig.5 Molecule collision rate in different number of simulation molecules

3 結(jié) 論

本文選取空氣中化學(xué)反應(yīng)最常見的三種氮氧類分子，使用DSMC對三種氣體在絕熱熱浴條件下進行了離解反應(yīng)模擬。離解度模擬結(jié)果顯示，對含有化學(xué)反應(yīng)的流場，DSMC模擬隨著流場溫度和氣體分子離解度的升高，正確模擬這一化學(xué)過程所要求的模擬分子數(shù)同樣呈上升趨勢，在本文的計算條件下，以氮氣為例，從離解度5%到離解度95%得到正確結(jié)果的模擬分子數(shù)閾值分別為250、550、850和1900，最高離解度與最低離解度模擬分子數(shù)閾值相差10倍左右。這一規(guī)律與Bird對普通流場DSMC模擬建議的統(tǒng)一20～30個模擬分子數(shù)不同，在含有化學(xué)反應(yīng)的流場計算中，準(zhǔn)確模擬要求的模擬分子數(shù)閾值與溫度（離解度）呈正相關(guān)，超出閾值之后，計算開銷增大，但是并不會給計算精度帶來明顯的變化。這一結(jié)果表明，在不同的流場物理條件下，DSMC要求的模擬分子數(shù)具有特異性，不能一概而論。

由于文中的討論是在絕熱封閉條件下進行的，所以在模擬的過程中，不會有新分子的加入。從理論上講，在真實流動的條件下，由于流場中會有新分子的流入，所以需要的網(wǎng)格模擬分子數(shù)會比表4列出的模擬分子數(shù)下限要少。但是，對于不同溫度下模擬分子數(shù)需求的規(guī)律是一致的，所以，本文得到的結(jié)論對于含有化學(xué)非平衡現(xiàn)象的流場模擬具有一定的參考價值。

［1］WU Q F，CHEN W F，et al.Rarefied gas dynamics［M］.National University of Defense Technology Publishing，2004.（In chinese）吳其芬，陳偉芳，等.稀薄氣體動力學(xué)［M］.湖南：國防科技大學(xué)出版社，2004.

［2］ALEXANDER F J，GARCIA A L，ALDER B J.Cell size dependence of transport coefficients in stochastic particle algorithms［J］.Phys.Fluids，1998，10（6）：1540－1542.

［3］BIRD G A.The direct simulation Monte Carlo method：current status and perspectives［A］.in：M.Mareschal（Ed.），Microscopic Simulations of Complex Flows，Plenum［M］，New York，1990.

［4］TZENG Pei－Yuan，LIU Ming－Ho.Influence of number of simulated particles on DSMC modeling of micro－scale Rayleigh－Be′nard flows［J］.InternationalJournalofHeatandMassTransfer，2005，45：2841－2855.

［5］BIRD G A，GALLIS M A，et al.Accuracy and efficiency of the sophisticated direct simulation Monte Carlo algorithm for simulating noncontinuum gas flows［J］.PhysicsofFluids，2009，21：017103.

［6］BOYD I D.Rotational and vibrational nonequilibrium effects in rarefied hypersonic flows［J］.JournalofThermophysicsand HeatTransfer.1990，4：478－484.

［7］DAAN Frenkel，BEREND Smit.Understanding molecular simulation：from algorithms to applications［M］.Academic Press，New York，2010.

［8］WU Q F，CHEN W F，et al.The DSMC method for the thermochemical nonequilibrium flows of high temperature rarefied gas［M］.National University of Defense Technology Publishing，1999.（In chinese）吳其芬，陳偉芳.高溫稀薄氣體熱化學(xué)非平衡流動的DSMC方法［M］.湖南：國防科技大學(xué)出版，1999.

［9］GOLDSWORTHY M J，MACROSSAN M N，ABDEL－JAWAD M M.Nonequilibrium reaction rates in the macroscopic chemistry method for direct simulation Monte Carlo calculations［J］.PhysicsofFluids，2007，19：066101.

［10］BORGNAKKE C，LARSEN P S.Statistical collision model for Monte Carlo simulation of polyatomic gas mixture［J］.Journal ofComputationalPhysics，1975，18：405－420.

［11］BIRD G A.Molecular gas dynamics and the direct simulation of gas flows［M］.Clarendon Press，Oxford，1994.

［12］BOYD I.Nonequilibrium chemistry modeling in rarefied hypersonic flows［A］.AIAA Thermophysics Conference［C］，33rd，Norfolk，VA；United State；28June－1July 1999.

［13］HERZBERG G.Molecular spectra and molecular structure Vol.spectra of diatomic molecules［M］.Princeton，1950.

［14］MALVIN H Kalos，PAULA A Whitlock.Mont Carlo method second，revised and enlarged edition［M］.Wiley，New York，2008.

［15］MICHAEL A Gallis，JOHN K Harvey.The modeling of chemical reactions and thermochemical nonequilibrium in particle simulation computations［J］.PhysicsofFluids，1998，10（6）：1344－1358.

［16］PARKER G J.Rotational and vibrational relaxation in diatomic gases［J］.Phys.Fluids，1959，2（4）：449－62.