王 焱，李 響

（遼寧工程技術(shù)大學(xué)電氣控制與工程學(xué)院，遼寧 葫蘆島125105）

1 引言

煤矸石是一種多孔介質(zhì)，暴露于空氣后，其空隙表面將吸附空氣中的氧氣。由于氧氣在矸石山表面或通過縫隙和大裂縫滲入矸石山內(nèi)部吸附潛伏，使得煤矸石緩慢氧化并釋放大量的熱，從而造成熱量累積，熱量的累積加速氧化并自燃。本文就是使用蒙特卡羅模擬方法與一些環(huán)境參量相結(jié)合建立起安全監(jiān)測模型。

在安全監(jiān)測模型中，蒙特卡羅 MC（Monte Carlo）模擬方法根據(jù)空氣分子在矸石山中的輸運(yùn)過程對(duì)分子的運(yùn)動(dòng)軌跡進(jìn)行模擬?？諝夥肿釉陧肥街械倪\(yùn)動(dòng)過程，本身就具有隨機(jī)性質(zhì)，一個(gè)由源點(diǎn)發(fā)出的粒子，在它的運(yùn)動(dòng)方向上，在哪一點(diǎn)碰撞是偶然的，但是具有一定的幾率分布；分子之間發(fā)生的碰撞又有各種幾率不盡相同的類型；碰撞后的能量和方向，又遵從一定的概率分布。采用MC方法可以將整個(gè)難解問題轉(zhuǎn)化為大量簡單的隨機(jī)實(shí)驗(yàn)。結(jié)合分子運(yùn)動(dòng)論，找出分子在煤矸石山當(dāng)中做無規(guī)則運(yùn)動(dòng)時(shí)哪些因素對(duì)分子運(yùn)動(dòng)產(chǎn)生影響，使用MC方法模擬煤矸石山當(dāng)中空氣分子的運(yùn)動(dòng)狀況，在矸石山內(nèi)部放入探測器觀察和記錄分子的運(yùn)動(dòng)路徑，從而可以得出矸石山中容易吸附氧分子的點(diǎn)，吸附氧分子越多說明越容易產(chǎn)生氧化，進(jìn)而產(chǎn)生自燃和爆炸的現(xiàn)象，針對(duì)這些潛在的危險(xiǎn)點(diǎn)提出相應(yīng)的防護(hù)措施，對(duì)于促使煤炭企業(yè)意識(shí)到矸石山爆炸的危害并主動(dòng)預(yù)防矸石山爆炸具有重要的現(xiàn)實(shí)意義。該過程是研究煤自燃機(jī)理、認(rèn)識(shí)自燃規(guī)律、建立預(yù)測和預(yù)防技術(shù)的基礎(chǔ)。

2 蒙特卡羅模擬法及其模擬過程

2.1 蒙特卡羅模擬法

蒙特卡洛模擬實(shí)際上就是模擬相當(dāng)數(shù)量的粒子在介質(zhì)中運(yùn)動(dòng)的狀況，使粒子運(yùn)動(dòng)的統(tǒng)計(jì)規(guī)律得以重現(xiàn)。不過，這種模擬不是用實(shí)驗(yàn)方法，而是利用數(shù)值方法和技巧，即利用隨機(jī)數(shù)來實(shí)現(xiàn)的。MC方法是以概率和統(tǒng)計(jì)理論方法為基礎(chǔ)的一種計(jì)算方法，將所求解的問題同一定的概率模型相聯(lián)系，用計(jì)算機(jī)實(shí)現(xiàn)統(tǒng)計(jì)模擬或抽樣，以獲得問題的近似解。

采用MC方法可以將整個(gè)難解決問題轉(zhuǎn)化為大量簡單的隨機(jī)實(shí)驗(yàn)。由于蒙特卡洛模擬方法具有程序結(jié)構(gòu)簡單、模擬過程靈活、不受問題條件限制、適于求解多維問題等優(yōu)點(diǎn)，因而被廣泛應(yīng)用。

2.2 蒙特卡羅模擬方法的步驟

首先，隨機(jī)選取一種初始分布，計(jì)算其能量，設(shè)這個(gè)能量為E1，一般來說可以產(chǎn)生三種擾動(dòng)使系統(tǒng)分布狀態(tài)發(fā)生變化：（1）使一個(gè)粒子平動(dòng)；（2）插入一個(gè)粒子；（3）刪除一個(gè)粒子，把變化后的分布能量設(shè)為E2。如果E2＜E1，則新的分布在統(tǒng)計(jì)求和中更重要，需要保留在統(tǒng)計(jì)平均的求和號(hào)中。如果E1＜E2，新的分布在統(tǒng)計(jì)求和中的貢獻(xiàn)不如原來的分布，此時(shí)，計(jì)算新舊兩種分布的玻爾茲曼因子的比值：

當(dāng)E1＜E2時(shí)，r是小于1的正數(shù)。用偽隨機(jī)數(shù)發(fā)生器產(chǎn)生一個(gè)小于1的正數(shù)a。如果a＜r，新的分布保留在統(tǒng)計(jì)平均的求和號(hào)中。如果r＜a，丟棄新的分布，恢復(fù)原來的分布，如此循環(huán)，構(gòu)成一個(gè)馬爾科夫鏈。

3 氧氣分子運(yùn)動(dòng)過程的蒙特卡洛模擬

目前，MC算法借鑒中子傳輸理論，初步建立了光在生物組織中的傳輸模型，在其他的粒子輸運(yùn)問題中有著廣泛的應(yīng)用。但是，在模擬分子運(yùn)動(dòng)過程理論中的應(yīng)用較少，主要是由于分子運(yùn)動(dòng)過程本身的多樣性和復(fù)雜性以及現(xiàn)有理論工具不成熟。在MC已被廣泛地成功應(yīng)用于各種粒子的輸運(yùn)問題的同時(shí)，利用它來解決氣體分子運(yùn)動(dòng)問題就成了一個(gè)很自然的想法。

本文主要是建立氧分子在矸石山中運(yùn)動(dòng)過程的模型。氣體分子在矸石山中的傳輸過程非常復(fù)雜，難以獲得較好的解析解，但可以利用蒙特卡羅方法模擬氧氣分子在矸石山中被吸附或者反射的隨機(jī)過程，進(jìn)而通過隨機(jī)數(shù)抽樣和統(tǒng)計(jì)來獲得氧分子的傳輸過程。采用蒙特卡羅方法建立分子運(yùn)動(dòng)模型，首先做以下假設(shè)：（1）流動(dòng)為自由分子流，分子之間無碰撞；（2）運(yùn)動(dòng)過程中分子與矸石體表面只發(fā)生反射和吸附作用；（3）在入口截面上，進(jìn)入矸石山的分子在位置上是均勻分布的，在角度上符合余弦分布規(guī)律。在確定分子的運(yùn)動(dòng)模型后，便可以根據(jù)對(duì)分子流追蹤的三步流程（入射、傳輸和統(tǒng)計(jì)）來構(gòu)建蒙特卡羅隨機(jī)抽樣模型。

3.1 氧氣分子運(yùn)動(dòng)過程

本文借助MC算法解決粒子的輸運(yùn)問題來解決分子運(yùn)動(dòng)問題，因?yàn)榉肿舆\(yùn)動(dòng)過程可被視為一種隨機(jī)散射的粒子輸運(yùn)模型，分子在矸石山中的傳輸行為就類似于中子的輸運(yùn)問題。具體講，可以用一種粒子的傳輸過程來模擬分子在矸石山體內(nèi)的傳播。粒子的數(shù)密度等價(jià)為分子運(yùn)動(dòng)勢能。圖1給出了氧分子在矸石山內(nèi)部隨機(jī)行走的示意圖。

Figure 1 Scheme of molecule paths inside a trapezoid hillock圖1 分子在梯形矸石山內(nèi)行走的示意圖

假設(shè)分子與煤矸石表面相互作用的基本參數(shù)ua、us及S（θ）完全是唯象的，且原則上可由實(shí)驗(yàn)直接測定。這正好可用來構(gòu)造一個(gè)概率分布模型，并可容易地實(shí)現(xiàn)概率分布抽樣。

分子運(yùn)動(dòng)學(xué)中MC方法的模擬過程可簡單歸結(jié)為：通過跟蹤大量的分子，記錄各個(gè)分子在山體中的行跡，從統(tǒng)計(jì)學(xué)可獲得分子行跡的統(tǒng)計(jì)規(guī)律，計(jì)算出山體中的分子勢能或其它參量的空間時(shí)間分布，以及矸石表面的反射和吸收等物理參量。圖2給出了該模擬過程的流程圖。

Figure 2 Smiplified flow diagram for the Monte Carlo smiulation of photon propagation圖2 Monte Carlo模擬流程圖

跟蹤的主要步驟為：

（1）按照入射條件確定起始跟蹤點(diǎn)和初始化分子流（權(quán)重、方向和位置）；

（2）確定氧氣分子下一次碰撞位置；

（3）確定分子在該位置上由于吸附而導(dǎo)致的權(quán)重衰減，選取適當(dāng)?shù)姆瓷湎辔缓瘮?shù)確定分子由于反射而導(dǎo)致的方向改變；

（4）如果分子權(quán)重小于設(shè)定的閾值或溢出時(shí)，則結(jié)束對(duì)該分子的跟蹤，進(jìn)而開始對(duì)下個(gè)分子跟蹤；

（5）否則重復(fù)（2）、（3）步驟，直至所設(shè)定的分子數(shù)都已被跟蹤完畢。

3.2 蒙特卡羅模擬方法

3.2.1 理論模型及坐標(biāo)系的建立

假設(shè)氧氣分子從矸石山外部進(jìn)入矸石山內(nèi)部是垂直進(jìn)入的，所以該問題是柱面對(duì)稱的，除了在矸石山外部表面建立一個(gè)直角坐標(biāo)系外，還建立一個(gè)柱面坐標(biāo)系來記錄模擬量，另外還建立一個(gè)移動(dòng)的球坐標(biāo)系來采樣分子流的傳播方向，如圖3所示。

3.2.2 相互作用點(diǎn)坐標(biāo)關(guān)系

設(shè)入射點(diǎn)為坐標(biāo)原點(diǎn)O（0，0，0），分子最初傳輸方向［7］為Z軸方向，用方向余弦表示為（0，0，1），權(quán)重因子W ＝1。分子運(yùn)動(dòng)的步長由平均自由程的累積概率分布決定。這里的反射點(diǎn)是隨機(jī)分布的，則平均自由程為u－1t（ut＝us＋ua，ut為總相互作用系數(shù)，us為反射系數(shù)，ua為吸附系數(shù)），所以平均自由程的累積分布可表示為：

其中，ζ是（0，1）上均勻分布的隨機(jī)數(shù)。

步長確定后，即可計(jì)算并確定分子流新位置的坐標(biāo)。如圖3所示，分子流在運(yùn)動(dòng)時(shí)，相鄰兩作用點(diǎn) Pm（xm，ym，zm）和 Pm＋1（xm＋1，ym＋1，zm＋1）之 間的坐標(biāo)關(guān)系為：

Figure 3 Coordinate system of simulate the process of molecular motion圖3 模擬分子運(yùn)動(dòng)過程的坐標(biāo)系

其中，α1、β1、γ1是Pm（xm，ym，zm）點(diǎn)分子流運(yùn)動(dòng)的方向余弦，lm是分子流在矸石山中運(yùn)動(dòng)到第m步時(shí)的步長，由式（1）確定。

3.2.3 化學(xué)吸附導(dǎo)致權(quán)重的衰減

分子流每運(yùn)動(dòng)一步，就應(yīng)計(jì)算在作用點(diǎn)處被矸石體表面吸附而導(dǎo)致的分子權(quán)重的衰減，同時(shí)在該點(diǎn)所處的柵元上記錄當(dāng)前權(quán)重因子被衰減的部分，得到該柵元上累計(jì)的分子的權(quán)重值。被矸石體吸收的權(quán)重ΔW 為：

修正后的權(quán)重為：

3.2.4 反射方向

權(quán)重衰減后的分子流在該作用點(diǎn)處還需經(jīng)歷一次反射事件。分子的反射由作用點(diǎn)處反射角的方向余弦描述。在球坐標(biāo)系下（圖3），分子經(jīng)過反射后，新的方向與反射前的方向夾角由方位角φ和反射角θ決定。根據(jù)Henyey－Greenstein相函數(shù)方程導(dǎo)出其抽樣值，再通過坐標(biāo)變換，分子流的運(yùn)動(dòng)新方向用方向余弦表示為：

從前面的討論中已看到，發(fā)生反射時(shí)，分子流的運(yùn)動(dòng)過程即自然終止。對(duì)于分子流仍在矸石山中傳輸?shù)錂?quán)重因子衰減到低于一定的閾值（如wt＝0.000 1）時(shí)，采用另一種降低方差的技巧，即所謂“俄羅斯輪盤”來對(duì)其進(jìn)行處理：給定一參數(shù)m（如m＝10），如隨機(jī)數(shù)ξ≤1／m，則分子流復(fù)活，其權(quán)重因子變?yōu)閙W，并繼續(xù)其運(yùn)動(dòng)過程；如隨機(jī)數(shù)ξ＞1／m，則權(quán)重因子為零，分子運(yùn)動(dòng)過程終止。

4 模擬及結(jié)果

為檢驗(yàn)上述數(shù)理模型和程序設(shè)計(jì)的正確性，選取一束分子包（其中含有10 000個(gè)氧氣分子）來對(duì)分子的運(yùn)動(dòng)軌跡進(jìn)行模擬，并對(duì)分子在不同的溫度和風(fēng)速等環(huán)境因素下的運(yùn)動(dòng)軌跡進(jìn)行比較。不同環(huán)境因素下的參數(shù)如表1所示。

Table 1 Parameters under different environmental factors表1 不同環(huán)境因素下的參數(shù)

進(jìn)行蒙特卡羅模擬時(shí)，隨機(jī)數(shù)產(chǎn)生的質(zhì)量會(huì)影響模擬的精度，為得到誤差較小的統(tǒng)計(jì)量，光子數(shù)至少為104個(gè)。以入射點(diǎn)為中心，計(jì)算不同半徑位置的分子密度流，在此模擬過程中可將分子束運(yùn)動(dòng)過程近似看成是光子流的運(yùn)動(dòng)過程，其密度流公式為：

其中 ，F(xiàn)是分子密度流，C是分子濃度，Ua是單元總體積，N是總光子數(shù)，V是單元格體積。

圖4為應(yīng)用表1所給的參數(shù)進(jìn)行模擬所得的在兩種環(huán)境因素影響下不同半徑處相對(duì)分子密度流。

Figure 4 Relative density of molecule in difference radius away from point of incidence圖4 距離入射點(diǎn)不同半徑處的相對(duì)分子密度流

可見，在球坐標(biāo)系中，在溫度和風(fēng)速的不同環(huán)境因素影響下，相對(duì)分子密度流隨著半徑的增大而減小。表1和圖4的結(jié)果表明，在不同環(huán)境下的吸收系數(shù)和散射系數(shù)越小，矸石山表面不同半徑位置處觀測到的分子密度流越大。

圖4中黑色點(diǎn)代表濕度影響下產(chǎn)生的密度流坐標(biāo)，～ 點(diǎn)代表在溫度影響下產(chǎn)生的密度流坐標(biāo)，×點(diǎn)代表風(fēng)速影響下的密度流坐標(biāo)。

5 結(jié)束語

本文初步應(yīng)用蒙特卡羅模擬方法對(duì)矸石山內(nèi)的分子運(yùn)動(dòng)過程建立基本模型，并對(duì)分子束在溫度、風(fēng)速、濕度等不同環(huán)境因素的影響下在矸石山內(nèi)部的傳播特性進(jìn)行研究。實(shí)驗(yàn)?zāi)M結(jié)果表明：

（1）相對(duì)的分子密度流隨著半徑的增大而減小。在距離入射中心較近處，相對(duì)分子密度流隨著半徑的增大下降很快。

（2）在濕度的影響下，吸收系數(shù)較大，它的相對(duì)分子密度流下降的趨勢相對(duì)比較快，而在風(fēng)速的影響下，吸收系數(shù)值較小，它的分子密度流下降的趨勢相對(duì)比較慢。

結(jié)果表明，矸石山內(nèi)部分子的密集度與不同環(huán)境因素下的吸收系數(shù)和散射系數(shù)都有密切聯(lián)系，所以，在對(duì)矸石山內(nèi)部分子運(yùn)動(dòng)建立基本模型的過程中，這些因素的影響是非常重要的。

［1］ Yuan Qin，Liu Yu－yan.Mathematical modeling and computer modeling process［J］.Journal of Chi zhou University，2010，24（6）：8－10.（in Chinese）

［2］ Deng Jun，Xu Jing－cai，Chen xiao－kun.Progress of research on the theory of coal spontaneous combustion mechanism and prediction［J］.Journal of Liaoning Technical University，2003，22（4）：455－145.（in Chinese）

［3］ Huang Wen－zhang.Study on mechanism and prevention technology of spontaneous combustion of coal gangue［D］.Chongqing：Chongqing University，2004.（in Chinese）

［4］ Li Ying.Monte Carlo particle transport in research on the key technology of pretreatment automatic modeling system calculation［D］.Hefei：Chinese Academy of Sciences，Hefei Institutes of Physical Science，2009.（in Chinese）

［5］ Xu Zhong－ji.monte carlo method［M］.Shanghai：Science and Technology Press，1980.

［6］ Zhang Lian－xing，Zheng Dian－chun，Chen Xue－feng.Monte Carlo simulation of particle dynamics in SF6－N2［J］.Journal of Harbin University of Science and Technology，2010，15（4）：89－93.

［7］ Peng Dong－qing，Li Hui，Xie Shu－sen.Monte Carlo method in tissue optics［Z］.Fuzhou：Laser Institute Fujian Normal University，2008.

［8］ Wang Ji－ren，Deng Cun－bao，Deng Han－zhong，et al.Microscopic mechanism of coal surface physical adsorption to oxygen molecules［J］.Coal Conversion，2007，30（4）：18－21.

［9］ Kernighan B W.Programming in C－A Tutorial［M］.Manh－attan：Murray Hill，1999.

［10］ Wu B T，Xiao D M.Predicting the dielectric strength of c－C4F8 200and SF6gas mixtures by monte carlo method［J］.Journal of Shanghai Jiaotong University，2007，（12E（1））：121－124.

［11］ Wang Li－h(huán)ong，Jacques S L.Monte Carlo modeling of light transport in multi－layered tissues standard C［J］.Computer Metheds and Programs in Biomedicine，1995，47（2）：131－146.

附中文參考文獻(xiàn)：

［1］ 袁琴，劉玉艷.數(shù)學(xué)建模與計(jì)算機(jī)建模過程［J］.池州學(xué)院學(xué)報(bào)，2010，24（6）：8－10.

［2］ 鄧軍，徐精彩，陳小坤.煤自燃機(jī)理及預(yù)測理論研究進(jìn)展［J］.遼寧工程技術(shù)大學(xué)學(xué)報(bào)，2003，22（4）：455－459.

［3］ 黃文章.煤矸石山自燃發(fā)火機(jī)理及預(yù)防技術(shù)研究［D］.重慶：重慶大學(xué)，2004.

［4］ 李瑩.蒙特卡羅粒子輸運(yùn)計(jì)算自動(dòng)建模系統(tǒng)預(yù)處理關(guān)鍵技術(shù)研究［D］.合肥：中國科學(xué)院合肥物質(zhì)科學(xué)研究院，2009.

［5］ 徐鐘濟(jì).蒙特卡洛方法［M］.上海：科學(xué)技術(shù)出版社，1980.

［6］ 張連星，鄭殿春，陳雪鋒.SF6－N2中粒子動(dòng)力學(xué)特性的蒙特卡洛仿真［J］.哈爾濱理工大學(xué)學(xué)報(bào)，2010，15（4）：89－93.

［7］ 彭東青，李暉，謝樹森.組織光學(xué)中的 Monte Carlo方法［Z］.福州：福建師范大學(xué)激光研究所，2008.

［8］ 王繼仁，鄧存寶，鄧漢忠，等.煤表面對(duì)氧分子物理吸附的微觀機(jī)理［J］.煤炭轉(zhuǎn)化，2007，30（4）：18－21.