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

（上海交通大學 航空航天學院，上海 200240）

0 引 言

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

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

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

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

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

1 DSMC熱化學非平衡模擬方法

1.1 概述

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

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

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

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

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

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

1.2 DSMC模擬離解反應

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

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

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

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

2 計算結果及分析

2.1 DSMC程序驗證

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

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

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

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

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

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

圖2 標準算例結果Fig.2 DSMC results of the standard examples

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

2.2 模擬分子數對DSMC模擬結果的影響

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結 論

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

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

［1］WU Q F，CHEN W F，et al.Rarefied gas dynamics［M］.National University of Defense Technology Publishing，2004.（In chinese）吳其芬，陳偉芳，等.稀薄氣體動力學［M］.湖南：國防科技大學出版社，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）吳其芬，陳偉芳.高溫稀薄氣體熱化學非平衡流動的DSMC方法［M］.湖南：國防科技大學出版，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.