魏寧斐，馮 焱，孫雯君，成永軍，董 猛，宋 伊，吳成耀，劉國庭

（蘭州空間技術(shù)物理研究所 真空技術(shù)與物理重點(diǎn)實(shí)驗(yàn)室，蘭州 730000）

0 引言

在真空技術(shù)中，混合氣體廣泛應(yīng)用于生物醫(yī)學(xué)[1]、航天器微型推力器[2]、分壓力質(zhì)譜計的校準(zhǔn)等。但是對于混合氣體條件下氣體分子之間的碰撞過程研究較少（平均自由程、碰撞頻率等）[3]，而研究混合氣體分子平均自由程在工程上有很重要的意義，如通過對克努曾數(shù)的計算從而正確判斷流動狀態(tài)等。相比單一氣體，混合氣體分子之間的碰撞除了考慮氣體的稀薄程度、壓力和溫度等參數(shù)外，還要考慮構(gòu)成混合氣體的組分以及各組分的濃度、混合物中每種組分的平均分子速度隨混合物的平均分子速度以及宏觀速度的不斷變化等[4]。所以對于混合氣體條件下分子之間碰撞過程的模擬要比對單一氣體的模擬更復(fù)雜。

目前大多數(shù)模擬研究是針對單一氣體分子內(nèi)部相互碰撞過程的，運(yùn)用的方法主要是分子氣體動力學(xué)。Einwohner等[5]、Dongari等[6]使用硬球模型對氣體分子的平均自由程和碰撞次數(shù)等微觀量進(jìn)行了模擬研究。李志剛等[7]研究了氣體分子與納米粒子的碰撞過程，進(jìn)一步研究了漫反射在氣體分子與納米粒子碰撞中產(chǎn)生的原因。何格等[8]利用Matlab對理想氣體下氣體的平均自由程等進(jìn)行了模擬計算，并驗(yàn)證了其算法中速度分布遵從麥克斯韋速度分布。李海濤等[9]基于分子氣體動力學(xué)方法，利用C語言自編模擬程序，對空氣的平均自由程、自由程分布等進(jìn)行了模擬研究。相比分子動力學(xué)模擬方法，直接模擬蒙特卡羅方法也常被人們用來研究分子的運(yùn)動、碰撞以及與壁面的相互作用[10-12]。直接模擬蒙特卡羅法首先由澳大利亞的Bird[10]提出，其基本原理是用大量的模擬分子模擬真實(shí)的氣體分子，每個模擬分子通常代表相當(dāng)數(shù)量的真實(shí)分子，須對分子的碰撞與運(yùn)動進(jìn)行解耦。White等[13]也運(yùn)用直接模擬蒙特卡羅方法對O2在單位時間單位體積內(nèi)的碰撞次數(shù)進(jìn)行了模擬計算，并與理論計算進(jìn)行比較，其結(jié)果一致性較好。直接模擬蒙特卡羅方法還被廣泛應(yīng)用于低地球軌道和深空中的衛(wèi)星和航天器的模擬、返回艙進(jìn)入大氣的模擬、太空環(huán)境下噴嘴和噴頭射流的模擬、微尺度流動的模擬和快速非平衡的氣體流動的模擬（激光燒蝕、蒸發(fā)、沉積）及全球大氣評估等方面[14-17]。

但是，這些模擬工作普遍存在兩方面的問題，其一，多數(shù)研究的是單一氣體分子的平均自由程，且對于氣體分子之間的碰撞處理均選用傳統(tǒng)的硬球模型（HS），即未考慮氣體分子的直徑是隨著其相對速度不斷變化這一現(xiàn)實(shí)因素；其二，多數(shù)是對單一溫度下碰撞過程進(jìn)行研究，未考慮寬溫度區(qū)間不同模型模擬結(jié)果的差異性。針對上述不足，本文采用直接模擬蒙特卡羅方法，選用可變硬球模型（VHS）以及可變軟球模型（VSS），分別對單一及二元混合氣體在平衡態(tài)下分子間的碰撞過程在溫度為300～20 000 K間進(jìn)行模擬研究，通過理論計算驗(yàn)證模擬結(jié)果的準(zhǔn)確性。同時將兩種模型下的模擬結(jié)果與HS的計算結(jié)果進(jìn)行比較。

1 物理模型

1.1 模擬空間模型的建立與邊界條件設(shè)置

為了模擬空間中單一氣體及混合氣體在平衡狀態(tài)下內(nèi)部分子的運(yùn)動與碰撞的微觀過程，包括分子單位時間單位體積內(nèi)的碰撞次數(shù)、平均自由程等，本次模擬建立了一個邊長為0.027 m，體積為1.96×10-5m3的正六面體計算區(qū)域，整個計算區(qū)域?yàn)橐粋€網(wǎng)格。本文模擬的分子數(shù)密度對于單一氣體及二元混合氣體均為1×1019個/m3，其中每個分子代表真實(shí)分子個數(shù)為2×108個。模擬時間步長為10-7s，總共模擬1 000個時間步長，模擬時間為10-4s，這樣設(shè)置的依據(jù)是用硬球模型下平均自由程計算公式對模擬分子質(zhì)量最小的N2進(jìn)行計算，其平均自由程約為0.160 m，同時模擬時間步長也可用N2的平均自由程與其最大平均速度（20 000 K）之比即1.2×10-4s，該值小于平均碰撞時間。另外采用無時間計數(shù)器（NTC）碰撞處理方法，可完全避免采用傳統(tǒng)時間計數(shù)（TC）方法時出現(xiàn)的問題，即在特別強(qiáng)激波前緣等高度非平衡情況下，概率與時間上間隔的關(guān)聯(lián)帶來的問題，偶然接收的小概率的碰撞可能導(dǎo)致網(wǎng)格時間往前推進(jìn)很長一段，以至于很多時間步長內(nèi)無碰撞計算。

模型初始化時，在模擬空間中產(chǎn)生一定數(shù)量的模擬分子，每個模擬分子代表一定數(shù)量的真實(shí)分子。直接模擬蒙特卡羅方法首先以隨機(jī)數(shù)的方式在模擬空間里生成每個模擬分子的初始位置坐標(biāo)，遵從麥克斯韋速度分布的初始速度，如式（1）所示：

式中：v為氣體分子運(yùn)動速度，m/s；k為玻耳茲曼常數(shù)，J/K；m為氣體分子質(zhì)量，kg；T為模擬溫度，K。

初始化溫度為300～20 000 K。圖1所示為本文模擬的物理模型，即單一或者混合氣體初始化后氣體分子空間分布圖。

圖1 模擬空間模型中的氣體分子三維空間分布圖Fig.1 Three-dimensional spatial distribution map of gas molecules in simulated space model

壁面邊界條件采用周期性邊界條件，這樣設(shè)置可以消除由于粒子進(jìn)入和離開表面所引入的表面效應(yīng)和邊界效應(yīng)。本文對正六面體計算域的壁面溫度均設(shè)置為300 K。相比于傳統(tǒng)的氣體分子動力學(xué)模擬方法，其存儲和參與運(yùn)算的數(shù)據(jù)量相對比較小，在相同的模擬空間模擬相同的分子數(shù)時，運(yùn)算時間大幅減小，模擬更接近于平時所遇到的情況，即分子數(shù)密度較大，模擬結(jié)果的準(zhǔn)確性高，重復(fù)性好。

1.2 二元碰撞模型

模擬過程中，氣體分子之間的碰撞采用VHS[18]。這是因?yàn)閭鹘y(tǒng)硬球分子具有固定不變的截面，而真實(shí)分子的碰撞截面是隨著相對速度的不同而改變的（隨著相對速度cr的增加而減小，同時隨著相對平動能εt的增大而減?。?。這是模擬真實(shí)氣體流動過程中必須要滿足的條件。這樣能保證氣體黏性系數(shù)μ與溫度T的關(guān)系和實(shí)際氣體一致，而不是傳統(tǒng)硬球模型下黏性系數(shù)μ與溫度T所呈現(xiàn)的線性關(guān)系。VHS模型假設(shè)分子具有硬球一樣的均勻散射概率，這使得計算碰撞后速度的取樣十分簡單。分析和數(shù)值研究的實(shí)際應(yīng)用表明，分子的散射概率可觀測到的變化很小[18]。

令碰撞截面正比于相對速度的逆冪次。引入總碰撞截面σT、分子直徑d、相對速度cr和相對平動能量εt=mrcr2/2 的參考值：σT，ref、dref、cr，ref和εt，ref，其中σT，ref、dref是相對速度為cr，ref時的值，可以定義VHS模型為：

其中：ξ是碰撞截面σT依賴于相對速度cr、平動能量εt的規(guī)律中的方冪，其值為常數(shù)。并且VHS模型散射角與硬球模型的規(guī)律相同。

同時也利用VSS模型做了對照模擬，VSS模型與VHS模型的不同之處在于考慮了動量截面與黏性截面的比。其直徑變化方式與VHS模型相同，但其散射角為：

式中：b為瞄準(zhǔn)距離，即質(zhì)心坐標(biāo)系中，尚未受到兩分子相互作用力影響的分子軌道間的最接近距離；d為氣體分子直徑；α為散射系數(shù)，其值一般在1～2之間。

1.3 不同二元碰撞模型下單一及混合氣體相關(guān)量計算

氣體分子自由程是克勞修斯于1859年提出的，其物理意義是一個分子相繼兩次碰撞之間所通過的路程。對于單一氣體，其分子平均自由程可以運(yùn)用基于分子速度分布規(guī)律推導(dǎo)的式（4）方便地計算：

式中：λ為分子的平均自由程，m；為分子平均熱運(yùn)動速率，m/s；f為碰撞頻率，次/s；n為分子數(shù)密度，m-3；d為氣體分子直徑，m。

當(dāng)擴(kuò)展到混合氣體時，由于其各組分的摩爾質(zhì)量不同導(dǎo)致了熱運(yùn)動速度的多樣性和復(fù)雜性，即每種組分不但有相對于自己組分的速度還有相對于其他組分的速度。這也給混合氣體的平均自由程計算帶來了挑戰(zhàn)。因?yàn)榛旌蠚怏w分子平均自由程是由分子之間的相對速度、平均碰撞截面、熱運(yùn)動速度等決定的。假設(shè)混合氣體由p組分和q組分構(gòu)成，可選取p組分分子作為測試分子，將q組分分子視為場分子。因此混合氣體中p、q兩組分分子之間的碰撞頻率為：

p組分分子的碰撞頻率可以通過對所有組分的碰撞求和得到，即：

所以p組分分子的平均自由程等于該組分的平均熱運(yùn)動速率除以碰撞頻率：

混合氣體分子平均自由程則可寫為：

方程式（5）～（8）只是混合氣體碰撞頻率和平均自由程的一般結(jié)果，即將分子作為傳統(tǒng)硬球模型來考慮。對于可變硬球模型（VHS），考慮s種VHS分子組成的混合氣體，則當(dāng)混合氣體處于平衡狀態(tài)時，則其中p分子與所有q分子之間的碰撞頻率可以寫成式（9）：

式中：dref為氣體分子參考直徑，m；T為溫度，K；mr為折合質(zhì)量，kg；ωpq為p、q分子的黏性系數(shù)；Tref為參考溫度，K。

由式（6）可得組分p的碰撞頻率：

則混合氣體組分p的平均自由程等于其平均熱速率除以其碰撞頻率，可得：

最后，混合氣體的全局平均自由程為：

式中：mp為p分子質(zhì)量，kg；mq為q分子質(zhì)量，kg；

因此，對于混合氣體來說，要想準(zhǔn)確獲取克努曾數(shù)進(jìn)而計算出其流態(tài)較為困難。更進(jìn)一步分析所要用到的一個結(jié)果是當(dāng)混合氣體處于平衡狀態(tài)時，平衡氣體中p組分和q組分之間單位時間單位體積內(nèi)的碰撞次數(shù)，即：

當(dāng)式（12）、式（13）中的ωpq=1/2時，可用作計算HS下混合氣體分子平均自由程及p組分和q組分之間單位時間單位體積內(nèi)的碰撞次數(shù)。而對于VSS來說，其相關(guān)量的計算方法和VHS計算方法相同，只是其在273 K時氣體分子的參考直徑不同。

2 模擬流程

本文依據(jù)Bird提出的直接模擬蒙特卡羅方法，運(yùn)用開源軟件OpenFOAM-2.4.0-MNF中的dsmc?Foam+對單一及二元混合氣體在平衡狀態(tài)下的微觀量進(jìn)行了模擬。選取模擬時間作為總體循環(huán)控制變量，模擬時間步長一般小于分子之間的平均碰撞時間，在每一個時間步長?t內(nèi)，首先認(rèn)為分子按照平衡狀態(tài)時麥克斯韋速度分布做勻速直線運(yùn)動，然后判斷其是否與壁面發(fā)生碰撞，若與壁面發(fā)生碰撞，則按照壁面邊界條件設(shè)置，本文壁面邊界條件為周期性邊界條件，即粒子到達(dá)壁面后消失，然后按照原速度從相對的壁面中發(fā)射出來。當(dāng)分子運(yùn)動到新的位置時，更新分子位置坐標(biāo)并記錄其速度大小、方向等量，按照一定的順序?qū)δM分子重新編號。然后，在每個網(wǎng)格中按照接收-拒絕原則判斷分子間是否發(fā)生碰撞，若分子之間發(fā)生碰撞，按照本文所選擇的VHS模型或者VSS模型計算碰撞后的速度和內(nèi)能等。當(dāng)模擬時間達(dá)到設(shè)定值后，跳出循環(huán)，自動結(jié)束本次模擬并對流場性質(zhì)進(jìn)行抽樣得出總的碰撞次數(shù)、平均自由程等值，同時也可以對模擬空間內(nèi)的壓力進(jìn)行統(tǒng)計。具體流程如圖2所示。

圖2 直接模擬蒙特卡羅方法流程圖Fig.2 Direct Simulation Monte Carlo method flow char

3 結(jié)果與分析

建立了長寬高均為為0.027 m，體積為1.96×10-5m3的正六面體模擬空間模型如圖1所示，并對N2、O2、Ar等單一氣體，以及N2∶O2、Ar∶O2、Ar∶N2分別為1∶1的混合氣體的分子平均自由程、單位時間單位體積內(nèi)的碰撞次數(shù)等在模擬溫度為300～20 000 K之間進(jìn)行了模擬計算（具體參數(shù)如表1和表2所列，其中表2中的數(shù)值可由參考文獻(xiàn)[18]查得），并將模擬計算結(jié)果與上述VHS模型理論計算結(jié)果相比較，驗(yàn)證模擬方法和模擬結(jié)果的準(zhǔn)確性。同時設(shè)置了在VSS下，單一氣體N2和N2∶Ar為1∶1時的平均自由程、單位時間單位體積內(nèi)的碰撞次數(shù)隨著溫度變化的對照模擬。并將上述兩種模型的模擬結(jié)果與傳統(tǒng)硬球模型的結(jié)果進(jìn)行對比，比較三種模型下結(jié)果的差異。

表1 單一及混合氣體分子數(shù)密度及模擬時間Tab.1 Single and mixed gas number density and simulation time

表2 不同模型下模擬分子特性Tab.2 Simulated molecular properties at different models

3.1 氣體分子自由程測試

統(tǒng)計了單一氣體以及二元混合氣體在不同溫度下模擬空間內(nèi)氣體分子的平均自由程，并與VHS模型下混合氣體全局平均自由程計算式（12）的理論計算結(jié)果進(jìn)行比較，如表3所列，從表中可以看出，隨著溫度的升高，單一氣體及二元混合氣體的平均自由程均增大，這是由于溫度升高其熱運(yùn)動速度增大從而使其碰撞截面減小所致。對于單一氣體，其相對偏差很小，最大只有0.08%，且隨著模擬溫度的升高沒有明顯變化。對于二元混合氣體，其模擬值和計算值之間的相對偏差與混合氣體兩組分相對分子質(zhì)量差值有關(guān)，當(dāng)兩組分的相對分子質(zhì)量差越小時，其模擬值與計算值越吻合。但是當(dāng)兩組分相對分子質(zhì)量相差較大時，如Ar∶N2=1∶1時模擬值與計算值之間的相對偏差較大，引起這種現(xiàn)象的原因是輕分子的碰撞具有較大的相對速率，而且這個相對速度要用在(σTcr)max中，于是較輕的分子接受率要高于較重的分子。

表3 VHS模型下單一及二元混合氣體分子平均自由程理論計算與模擬值比對表Tab.3 Simulation and calculated mean free path of pure and binary gas mixture at VHS model

本文還利用VSS模型對N2及N2與Ar的混合氣體的平均自由程進(jìn)行了模擬計算，并利用式（12）進(jìn)行了理論計算，詳細(xì)計算結(jié)果與模擬結(jié)果的對比如表4所列。并與用HS計算的結(jié)果進(jìn)行對比，如圖3和圖4所示，發(fā)現(xiàn)在HS下分子數(shù)密度相同時單一及混合氣體分子平均自由程在寬溫度區(qū)間有很大的誤差，這是由于HS的分子直徑是一個定值，不隨溫度變化。而對于VHS和VSS而言，考慮了分子直徑隨著溫度升高會變小這一現(xiàn)實(shí)情況，其平均自由程隨著模擬溫度的升高而增大。但是可變軟球模型的值要更大，是因?yàn)槠淇紤]了動量截面與黏性截面之比，從文獻(xiàn)[19]中可知，N2、O2在VSS模型下黏性系數(shù)值在溫度為300～15 000 K間與實(shí)驗(yàn)數(shù)值符合很好，所以VSS模型與真實(shí)氣體條件更接近。因此，本文中VSS模型下平均自由程模擬結(jié)果應(yīng)更符合實(shí)際情況。

圖3 N2在不同溫度不同模型下的氣體分子平均自由程Fig.3 The mean free path of N2under different temperatures and different models

圖4 N2∶Ar=1∶1在不同溫度不同模型下的氣體分子平均自由程Fig.4 The mean free path of N2∶Ar=1∶1under different temperatures and different models

表4 VSS模型下單一及二元混合氣體分子平均自由程理論計算與模擬值比對表Tab.4 Simulation and calculated mean free path of pure and binary gas mixture at VSS model

3.2 氣體分子碰撞次數(shù)測試

表5統(tǒng)計了三種單一氣體以及三種混合氣體在總的分子數(shù)密度不變時VHS模型下單位時間單位體積內(nèi)的碰撞次數(shù)隨模擬溫度的變化關(guān)系。表6對單一氣體N2和N2∶Ar=1∶1的混合氣體在VSS模型下進(jìn)行了對照模擬。

表5 VHS模型單一氣體及二元混合氣體單位時間單位體積內(nèi)分子的碰撞次數(shù)Tab.5 Simulated collision rates compared with analytical values for pure and binary gas mixtures at VHS

表6 VSS模型單一氣體及二元混合氣體單位時間單位體積內(nèi)分子的碰撞次數(shù)Tab.6 Simulated collision rates compared with analytical values for pure and binary gas mixtures at VSS

運(yùn)用式（13）分別對VHS和VSS兩種模型下單位時間單位體積內(nèi)的碰撞次數(shù)進(jìn)行了計算，并對模擬值和理論計算值進(jìn)行了比較，由表5及表6可看出，理論值與模擬值相差不大，在概率波動范圍內(nèi)。在同一模型下其單位時間單位體積內(nèi)總的碰撞次數(shù)與單一氣體的分子質(zhì)量或者與混合氣體的折合質(zhì)量有關(guān)，即在相同的分子數(shù)密度、相同的溫度條件下，其單位時間單位體積內(nèi)碰撞次數(shù)隨著分子質(zhì)量或者折合質(zhì)量的減小而增大。在相同的分子數(shù)密度下，同種氣體分子之間的碰撞次數(shù)隨著溫度的升高也不斷增大，這種現(xiàn)象符合布朗運(yùn)動的特征。

將傳統(tǒng)HS的計算結(jié)果（其計算方法是令式（13）中的ωpq=1/2，且dN2=3.754×10-10m，dAr=3.628×10-10m）與VHS和VSS兩種模型下的模擬結(jié)果進(jìn)行比較，具體如圖5圖6所示。對于單一氣體及混合氣體，其HS的結(jié)果與溫度的二分之一次方成正比，當(dāng)模擬溫度大于300 K時，其值大于VHS和VSS的結(jié)果。其大小從VSS到VHS再到HS由小變大，這說明VHS模型考慮了真實(shí)的分子碰撞截面是隨著相對速度的不同而改變，黏性系數(shù)隨溫度不同也接近于真實(shí)氣體，VSS模型除了VHS所考慮的因素外還考慮了動量截面與黏性截面的比，使得其結(jié)果更接近于真實(shí)條件，所以在混合氣體條件下，利用VSS更加準(zhǔn)確[18]。

圖5 N2在不同模型下單位時間單位體積內(nèi)分子碰撞次數(shù)Fig.5 Collision rates of N2at different models

圖6 Ar∶N2=1∶1在不同模型下單位時間單位體積內(nèi)分子碰撞次數(shù)Fig.6 Collision rates of Ar∶N2=1∶1 at different models

從模擬結(jié)果與計算結(jié)果之間的相對偏差可知，針對單一氣體及二元混合氣體在平衡狀態(tài)下氣體分子間碰撞過程的直接模擬蒙特卡羅方法是準(zhǔn)確的。相比于利用分子動力學(xué)模擬計算平均自由程的方法，直接模擬蒙特卡羅方法所用的模擬時間更短、單一氣體的模擬計算精度更高，李海濤等[9]對空氣模擬計算的相對誤差最大為12.1%。能夠用來計算混合氣體的全局平均自由程，只要混合氣體組分間相對分子質(zhì)量相差較小，其相對偏差不超過0.8%。相比于何格等利用傳統(tǒng)HS對He在298 K時平均自由程的模擬，本文對分子之間的碰撞處理運(yùn)用了VHS及VSS，模擬溫度為300～20 000 K。同時對二元混合氣體情況下的分子平均自由程、單位時間單位體積內(nèi)的碰撞次數(shù)也進(jìn)行了模擬計算，因此對處理兩種及兩種以上混合氣體全局平均自由程具有指導(dǎo)及借鑒意義。

4 結(jié)論

對于單一氣體，氣體分子平均自由程及碰撞次數(shù)的模擬結(jié)果與理論計算結(jié)果的相對偏差很小，對于相對分子質(zhì)量接近的二元混合氣體，其相對偏差亦很小。隨著兩組分相對分子質(zhì)量之差變大，模擬與計算值的相對偏差也變大，驗(yàn)證了VHS及VSS、壁面周期性邊界條件設(shè)置、模擬分子數(shù)密度1×1019m3等條件下利用直接模擬蒙特卡羅方法的可行性與準(zhǔn)確性。將VHS與VSS的模擬結(jié)果與HS的計算結(jié)果進(jìn)行比較，發(fā)現(xiàn)HS下的氣體分子的平均自由程在寬溫度區(qū)間有較大誤差，氣體分子碰撞次數(shù)在模擬溫度大于300 K時也遠(yuǎn)高于其他兩種模型的結(jié)果。所以在處理較高溫度時氣體分子之間的碰撞過程選用VHS或者VSS，結(jié)果會更準(zhǔn)確。在分子水平研究了單一及二元混合氣體的碰撞過程，為后續(xù)二元混合氣體在總分子數(shù)密度相同下，按照不同比例混合時氣體的平均自由程、單位時間單位體積內(nèi)的碰撞次數(shù)模擬研究，以及更多組分的混合氣體模擬研究提供基礎(chǔ)。有助于解決一些實(shí)際工程問題，例如通過計算混合氣體全局平均自由程從而對混合氣體克努曾數(shù)進(jìn)行計算，進(jìn)而能準(zhǔn)確判斷氣體的流動狀態(tài)。