煤炭自然發(fā)火介尺度分析:從表征體元宏觀模型到孔隙微觀模型

梁運(yùn)濤,王樹剛,蔣 爽,胡沛裕,林 琦,宋雙林,

(1.煤科集團(tuán)沈陽研究院有限公司 煤礦安全技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 撫順 113122; 2.大連理工大學(xué) 建設(shè)工程學(xué)部，遼寧 大連 116024; 3.大連民族大學(xué) 土木工程學(xué)院，遼寧 大連 116600)

煤炭作為我國(guó)主體能源，其安全生產(chǎn)直接關(guān)系到國(guó)民經(jīng)濟(jì)的平穩(wěn)發(fā)展[1]。煤礦火災(zāi)嚴(yán)重制約煤炭企業(yè)的安全生產(chǎn)，其中煤自燃是引發(fā)煤礦火災(zāi)的最主要原因[2-3]。煤炭自然發(fā)火研究一直是煤礦安全領(lǐng)域世界范圍的研究熱點(diǎn)與難點(diǎn)。

國(guó)內(nèi)外學(xué)者較早就從分子層面的物理化學(xué)過程開始研究和探索煤炭自然發(fā)火的原因。對(duì)于煤氧復(fù)合的微觀反應(yīng)歷程，通常采用吸附理論和氧化反應(yīng)來闡述，認(rèn)為煤自燃是煤氧之間吸附反應(yīng)和化合反應(yīng)同時(shí)作用的結(jié)果[4-10]。隨著現(xiàn)代科技的發(fā)展，紅外光譜分析技術(shù)、電子順磁共振光譜分析技術(shù)、核磁共振波譜分析技術(shù)等實(shí)驗(yàn)方法在煤自燃的研究中得到了日益廣泛的應(yīng)用，在分子層面對(duì)煤自燃過程開展了更深入的研究。基于分子結(jié)構(gòu)理論，煤中有機(jī)大分子側(cè)鏈基團(tuán)和低分子化合物等活性基團(tuán)是誘導(dǎo)煤自燃的物質(zhì)，活性基團(tuán)反應(yīng)特性可為闡釋煤炭微觀自燃機(jī)理和自燃過程提供科學(xué)依據(jù)[11-15]。

除分子層面的物理化學(xué)反應(yīng)影響外，煤巖自身的微觀物理結(jié)構(gòu)特征，如孔隙結(jié)構(gòu)、表觀結(jié)構(gòu)等，直接影響低溫氧化的發(fā)展過程[16-19]。基于顆粒聚團(tuán)尺度的研究發(fā)現(xiàn)，不同粒度的煤樣，比表面積不同，與煤樣發(fā)生氧化反應(yīng)的速度也不相同[20-22]，但是現(xiàn)有研究在關(guān)聯(lián)微觀特征與宏觀特性方面存在不足。

在煤自然發(fā)火宏觀模擬方面，通常采用多孔介質(zhì)理論進(jìn)行研究。以煤堆積體或煤自熱實(shí)驗(yàn)裝置為研究對(duì)象，考察介質(zhì)孔隙率、粒徑大小、外部風(fēng)速等諸多因素對(duì)煤堆自燃過程的影響，分析煤自熱過程中流場(chǎng)、溫度場(chǎng)的分布和熱點(diǎn)的遷移規(guī)律[23-28]。此外，還有針對(duì)采空區(qū)自然發(fā)火的多場(chǎng)耦合模擬研究[29-30]，巷道松散煤體自燃三維多場(chǎng)耦合模擬研究[31]，將流場(chǎng)、溫度場(chǎng)、濃度場(chǎng)、應(yīng)力場(chǎng)等進(jìn)行耦合，從礦井尺度上研究煤自燃高溫區(qū)域的發(fā)生及發(fā)展規(guī)律[32-34]。

綜上所述，在煤自燃領(lǐng)域，既存在著從分子層面活性基團(tuán)反應(yīng)著眼的機(jī)理性研究，也存在著從煤顆粒聚團(tuán)尺度的孔隙率、比表面積、粒徑大小等角度出發(fā)的煤自燃氧化反應(yīng)速度的研究，還存在著煤堆積尺度、采空區(qū)尺度、礦井尺度等的煤自燃火災(zāi)發(fā)生蔓延規(guī)律的研究，即煤炭自然發(fā)火存在著明顯的多尺度現(xiàn)象。

煤炭自然發(fā)火源于煤分子尺度下的物理吸附、化學(xué)吸附和化學(xué)反應(yīng)，但是否形成非控制燃燒主要取決于堆積煤體尺度下的高溫區(qū)發(fā)生發(fā)展特性。因此，采用合適的方法將煤自然發(fā)火的微觀機(jī)理與宏觀變化特征相關(guān)聯(lián)是煤自然發(fā)火領(lǐng)域的重要研究方向。在前期研究中，筆者[35]借鑒復(fù)雜系統(tǒng)理論和多尺度科學(xué)研究思路，提出煤自燃介尺度方法，定義了煤炭自然發(fā)火研究領(lǐng)域的介尺度(圖1)，并探討了描述型、關(guān)聯(lián)型和極值型3種多尺度方法[36-37]在煤炭自然發(fā)火介尺度研究中的適用性。

筆者基于文獻(xiàn)[35]提出的煤自燃介尺度方法，進(jìn)一步闡述在介尺度II所在的堆積態(tài)煤體層次上所采用的建模方法。首先分析了尺度之間的關(guān)聯(lián)和信息傳遞，從數(shù)學(xué)建模角度明確介尺度II(圖1中的堆積態(tài)煤體層次)的本質(zhì)和作用。其次論述了表征體元尺度宏觀數(shù)學(xué)模型，結(jié)合煤自燃特性提煉出需要由孔隙尺度模型提供的重要參數(shù)。然后描述了孔隙尺度微觀數(shù)學(xué)模型，用實(shí)例介紹了表征體元的物理建模和提取方法。最后提出從孔隙尺度到表征體元尺度參數(shù)上聯(lián)的方法。

圖1 煤自燃領(lǐng)域內(nèi)的兩個(gè)介尺度[35]Fig.1 Two kinds of mesco-scales in coal spontaneous combustion[35]

1 堆積態(tài)煤體層次內(nèi)的信息傳遞

從反應(yīng)機(jī)理角度出發(fā)，介尺度II所在的堆積態(tài)煤體層次包含3個(gè)物理尺度，即煤顆粒尺度、顆粒聚團(tuán)尺度和煤堆(或采空區(qū)、礦井)尺度，如圖1所示。其中煤顆粒尺度和煤堆尺度稱為邊界尺度，特征易于表征和分析，研究相對(duì)成熟;煤顆粒聚團(tuán)尺度介于煤顆粒尺度和煤堆尺度之間，是連接2個(gè)邊界尺度的橋梁，稱為介尺度。

不同尺度之間相互關(guān)聯(lián)，存在信息傳遞。例如，煤顆粒聚團(tuán)內(nèi)部結(jié)構(gòu)特征決定了煤堆的局部孔隙率、滲透率等宏觀特征參數(shù)。同樣的，煤顆粒聚團(tuán)內(nèi)顆粒表面特征、氣固反應(yīng)引起的骨架變化等決定了煤堆內(nèi)部宏觀對(duì)流與擴(kuò)散的物理特征參數(shù)。由此可見，尺度之間正確的信息傳遞是多尺度模擬的前提，而且煤自燃領(lǐng)域尺度信息的傳遞主要是以“尺度上聯(lián)”形式實(shí)現(xiàn)，即煤堆尺度下的數(shù)學(xué)模型準(zhǔn)確性依賴于單煤顆粒尺度下信息的合理獲取。

顆粒聚團(tuán)的介尺度研究可以將單煤顆粒自燃特性與煤堆自然發(fā)火過程聯(lián)系在一起。物理結(jié)構(gòu)上煤堆由顆粒聚團(tuán)組成，顆粒聚團(tuán)則由單煤顆粒組成。從數(shù)學(xué)建模角度，為便于運(yùn)用數(shù)學(xué)分析的連續(xù)函數(shù)工具，宏觀上一般忽略復(fù)雜的煤堆積體內(nèi)部孔隙結(jié)構(gòu)，將煤堆積體視為滿足基于連續(xù)性介質(zhì)假定的多孔介質(zhì)，其最小單元或研究單元的質(zhì)點(diǎn)是由復(fù)雜的孔隙和骨架結(jié)構(gòu)組成的煤顆粒聚團(tuán)，即，顆粒聚團(tuán)是組成煤堆的最小單元，聚團(tuán)內(nèi)部結(jié)構(gòu)變化對(duì)物理過程的影響可以忽略，聚團(tuán)整體的宏觀物理量變化反映了聚團(tuán)內(nèi)眾多煤顆粒在自燃過程中相互影響的凈作用效果?；诖擞^點(diǎn)，可以認(rèn)為煤顆粒聚團(tuán)的本質(zhì)就是煤堆(介質(zhì))分析中的表征體元。

以表征體元為最小單位的宏觀模型稱為表征體元尺度宏觀模型。在表征體元尺度下，忽略聚團(tuán)內(nèi)結(jié)構(gòu)變化，所采用的數(shù)值模擬方法是一種基于宏觀數(shù)學(xué)模型的常規(guī)計(jì)算方法。表征體元尺度是一種比孔隙尺度大得多的尺度，在表征體元尺度下多孔介質(zhì)內(nèi)宏觀物理量不再受到微觀孔隙結(jié)構(gòu)的影響，可以被視為一種連續(xù)介質(zhì)。以孔隙率的確定為例，從多孔介質(zhì)物理結(jié)構(gòu)中取出極小一塊區(qū)域測(cè)定其孔隙率，由于孔隙結(jié)構(gòu)的隨機(jī)性與特異性，當(dāng)測(cè)定區(qū)域變化時(shí)，孔隙率數(shù)值將出現(xiàn)震蕩。隨著測(cè)定區(qū)域體積增大，這種震蕩逐漸減小，直到某一體積，這種震蕩消失，孔隙率不再受到孔隙結(jié)構(gòu)隨機(jī)變化的影響，此體積即稱為表征體元。除孔隙率外，還可以定義出以表征體元為基本單位的宏觀或表觀流體物理量，如宏觀密度、速度、溫度和壓力等[38-39]。

相對(duì)于表征體元宏觀模型，以煤顆粒聚團(tuán)(表征體元)為研究對(duì)象所建立的物理和數(shù)學(xué)模型，本文之稱為孔隙尺度微觀模型。孔隙尺度微觀模擬考慮聚團(tuán)內(nèi)部孔隙結(jié)構(gòu)、單煤顆粒相互影響、氣固相互作用、煤顆粒熱解變化等過程，需要建立孔隙或者骨架結(jié)構(gòu)的物理模型。

基于上述分析可知，文獻(xiàn)[35]提出的煤顆粒聚團(tuán)這一介尺度做為橋梁將單煤顆粒和煤堆積體相關(guān)聯(lián)。煤顆粒聚團(tuán)既是宏觀尺度模型的最小單元或表征體元，又是微觀尺度模型的物理區(qū)域。出于計(jì)算效率和模型有效性的考慮，針對(duì)不同尺度建模往往選用不同的數(shù)學(xué)模型和計(jì)算方法[38]。

2 表征體元尺度宏觀數(shù)學(xué)模型

煤巖體是一種典型的多孔介質(zhì)，表征體元尺度宏觀控制方程包括連續(xù)性方程、動(dòng)量守恒方程、能量守恒方程和濃度方程。傳統(tǒng)的宏觀尺度數(shù)值計(jì)算方法往往基于連續(xù)性假設(shè)[40-41]，將煤巖體中的氣體視為不可壓縮流體，采用固定孔隙率、滲流達(dá)西定律、局部熱平衡等假設(shè)。然而，煤自然發(fā)火過程會(huì)導(dǎo)致煤巖體溫度不斷升高，一方面熱應(yīng)力會(huì)引起煤巖體結(jié)構(gòu)的變形，另一方面水分蒸發(fā)、原始吸附氣體解析、煤熱解等因素會(huì)導(dǎo)致煤的孔隙結(jié)構(gòu)更為發(fā)達(dá)[13,42-43]，從而引起局部孔隙率的變化。這種瞬時(shí)孔隙率可以根據(jù)實(shí)驗(yàn)數(shù)據(jù)擬合出公式供宏觀控制方程使用[31]。考慮這一特征的多孔介質(zhì)質(zhì)量守恒方程[44]為

(1)

其中，ε為孔隙率，與熱解溫度、空間位置有關(guān);ρ為密度，kg/m3;V為表觀速度，m/s;下標(biāo)f為流體。在滿足質(zhì)量守恒的同時(shí)，煤巖體內(nèi)的氣體滲流還遵循動(dòng)量守恒定律?？紤]到煤礦井下存在著諸如裂隙、采空區(qū)等大孔隙率的區(qū)域，動(dòng)量方程可采用Brinkman-Forchheimer 擴(kuò)展形式達(dá)西阻力模型，則三維非穩(wěn)態(tài)動(dòng)量方程[41]可表示為

(2)

其中，μ為流體動(dòng)力黏度，Pa·s;K為滲透率，m2;cF為慣性系數(shù);p為壓強(qiáng)，Pa。式(2)右邊第1項(xiàng)為壓力項(xiàng)，第2項(xiàng)為宏觀黏滯阻力項(xiàng)，第3項(xiàng)為微觀黏滯阻力(Darcy項(xiàng))與慣性力項(xiàng)。

基于能量守恒，煤巖體內(nèi)滲流過程中空氣與煤巖固體各自的能量方程可以寫為

ah(Ts-Tf)

3)

(1-ε)qs+F-ah(Ts-Tf)

(4)

其中，下標(biāo)s和f分別為固體和流體;T為溫度，K;C和Cp分別為固體比熱和流體定壓比熱，J/(kg·K);t為時(shí)間，s;h為對(duì)流換熱系數(shù)，W/(m2·K);a為流固相互作用面積密度，m-1;λ為導(dǎo)熱系數(shù)，W/(m·K);F為輻射通量向量，W/m2;qs為固體生熱項(xiàng)，W/m3。如果考慮煤巖體內(nèi)滲流速度較低，空氣與煤巖體可以進(jìn)行充分換熱，假定流體溫度與煤巖體溫度一致(即Tf=Ts=T)，則傳熱方程[41]可表示為

(5)

其中，下標(biāo)m為體積平均。式(5)左邊第1項(xiàng)為單位時(shí)間內(nèi)的能量積累效應(yīng)，第2項(xiàng)為對(duì)流傳熱項(xiàng);右邊第1項(xiàng)分別為導(dǎo)熱傳熱量，第2和3項(xiàng)為內(nèi)部產(chǎn)熱量和輻射熱量。qs與煤氧反應(yīng)過程有關(guān)，可計(jì)算[41]如下:

其中，c(O2)為氧氣物質(zhì)的量濃度，mol/m3;k0為指前因子，1/s;E為活化能,J/mol;R為氣體常數(shù)，J/(mol·K);n為反應(yīng)級(jí)數(shù);ΔH為煤氧化學(xué)反應(yīng)過程中生熱量，J/mol O2，可通過煤樣實(shí)驗(yàn)獲得;sv為形狀因子。當(dāng)煤巖體局部溫度高于一定數(shù)值后，輻射通量F不可忽略。輻射通量計(jì)算有各種不同的方法，對(duì)于煤巖體可以選用較為簡(jiǎn)單的近似處理方法，輻射傳熱被近似處理成為一個(gè)擴(kuò)散過程，類似于熱傳導(dǎo)，計(jì)算[45]如下:

F=βRT3T=λrT

式中，βR為常數(shù)，W/(m·K4);λr為輻射導(dǎo)熱系數(shù)，W/(m·K)。

煤巖體內(nèi)氣體物質(zhì)濃度方程[41]可表示為

式中，Yi為物質(zhì)濃度，kg/m3;Di為動(dòng)力彌散系數(shù)，m2/s;Si為生成/耗散項(xiàng)，kg/(m3·s)。

在上述方程中，孔隙率初始分布與空間位置有關(guān)，比如采空區(qū)采動(dòng)影響下，由于壓實(shí)程度不同而呈現(xiàn)隨空間變化的特征[33]。對(duì)于煤巖體的緩慢升溫過程，其局部孔隙率變化與煤種、熱解溫度有關(guān)，實(shí)踐中可采用實(shí)驗(yàn)的方法，在給定煤樣的情況下進(jìn)行升溫?zé)峤鈱?shí)驗(yàn)，通過電鏡掃描表面或者工業(yè)CT掃描內(nèi)部結(jié)構(gòu)的方法獲取局部孔隙率與熱解溫度的關(guān)系式。依次求解連續(xù)性方程、動(dòng)量方程和能量方程，當(dāng)溫度場(chǎng)發(fā)生變化時(shí)更新局部孔隙率，從而獲取孔隙率變化時(shí)的流場(chǎng)及溫度場(chǎng)。式(3)和(4)中的對(duì)流換熱系數(shù)h、式(2)中的滲透系數(shù)K和慣性系數(shù)cF可以通過傳統(tǒng)的實(shí)驗(yàn)測(cè)試/經(jīng)驗(yàn)公式獲取，也可以通過孔隙尺度的微觀方法模擬獲取。

3 孔隙尺度微觀模型

3.1 孔隙尺度數(shù)學(xué)模型

孔隙尺度下，煤顆粒聚團(tuán)內(nèi)部結(jié)構(gòu)具有強(qiáng)烈非均質(zhì)特性，基于連續(xù)性假定的數(shù)值方法不能有效地解決此問題，因此需要基于離散方法的數(shù)值算法，將多孔介質(zhì)表示為離散單元，比如孔隙網(wǎng)絡(luò)模型或格子Boltzmann方法(LBM)[38]。在煤微觀結(jié)構(gòu)的滲流研究中，由于多孔介質(zhì)的孔隙極度不規(guī)則，當(dāng)介質(zhì)骨架作為流場(chǎng)邊界時(shí)，傳統(tǒng)流體計(jì)算方法在模擬時(shí)對(duì)于邊界處理較為困難，而LBM的粒子背景使之可以方便、高效地處理這種具有復(fù)雜邊界的流場(chǎng)[39]，并得到了成功應(yīng)用[46]。孔隙尺度微觀數(shù)學(xué)模型包含用于計(jì)算孔隙內(nèi)流動(dòng)的標(biāo)準(zhǔn)LBM模型和計(jì)算含內(nèi)熱源的熱格子玻爾茲曼模型，具體方程[38]如下:

(9)

gi(x+eiδt,t+δt)-gi(x,t)=

(10)

格子Boltzmann方法的邊界處理需要根據(jù)已知的宏觀條件確定相應(yīng)邊界點(diǎn)上相應(yīng)分布函數(shù)的取值，目前已發(fā)展出啟發(fā)式格式、動(dòng)力學(xué)格式、外推格式等多種邊界處理格式，其中啟發(fā)式格式中的各類反彈格式操作簡(jiǎn)單，特別適用于復(fù)雜的不規(guī)則邊界的處理，使得LBM在模擬包含復(fù)雜孔隙結(jié)構(gòu)的多孔介質(zhì)流動(dòng)及傳熱方法具有較大優(yōu)勢(shì)并得到了大量的應(yīng)用[38-39,47]。

3.2 孔隙尺度物理建模

3.2.1孔隙結(jié)構(gòu)生成

孔隙尺度的模擬需要了解固體骨架的詳細(xì)信息，一般有兩種構(gòu)造煤巖體多孔介質(zhì)物理結(jié)構(gòu)的方法。一種是人工生成方法，有球體沉降法[48]、硬球Monte-Carlo方法[49]、分?jǐn)?shù)布朗運(yùn)動(dòng)法[50]、四參數(shù)隨機(jī)生長(zhǎng)(QSGS)方法[51]等，其中QSGS方法在孔隙尺度的LBM研究中應(yīng)用較多[46,52-55]。QSGS的基本方法是首先在空間中布置一定數(shù)量的生長(zhǎng)核，然后在生長(zhǎng)核的基礎(chǔ)上以不同的概率向空間不同方向生長(zhǎng)，重復(fù)此過程直到所在相的孔隙率滿足要求。圖2是利用該方法生成的二維多孔介質(zhì)結(jié)構(gòu)。

圖2 使用QSGS方法生成的多孔介質(zhì)結(jié)構(gòu)(200像素×200像素)[56] Fig.2 Schematics of the generated porous structures using QSGS method in 200 grids×200 grids[56]

構(gòu)造煤巖體多孔介質(zhì)物理結(jié)構(gòu)的另一種方法是使用CT顯微斷層掃描技術(shù)對(duì)其進(jìn)行掃描與重構(gòu)。微米量級(jí)分辨率的工業(yè)CT掃描技術(shù)在巖石孔隙成像領(lǐng)域有廣泛應(yīng)用，但在煤自然發(fā)火領(lǐng)域應(yīng)用較少。一種立式旋轉(zhuǎn)的工業(yè)CT設(shè)備基本原理如圖3所示[57]，樣品被固定于CT工作倉內(nèi)的可旋轉(zhuǎn)操作臺(tái)上，微焦點(diǎn)X光源固定，樣品沿水平面方向勻速旋轉(zhuǎn)360°以獲得整體樣品的透視圖像，掃描完成后，對(duì)圖像進(jìn)行分解獲得所有斷層的二維層析圖像。本文以活性炭試樣的掃描圖像為例，說明根據(jù)CT顯微斷層掃面圖像計(jì)算多孔介質(zhì)孔隙率及提取表征體元的基本步驟。

圖3 CT掃描成像原理示意[57]Fig.3 Schematic of scanning and imaging[57]

利用工業(yè)CT對(duì)多孔介質(zhì)掃描后可得到沿某一方向有序排列的顯微斷層圖像。這些圖像是根據(jù)多孔介質(zhì)內(nèi)孔隙與固體骨架對(duì)X射線吸收率不同得到的0～255階像素的灰度圖，其中圖4(a)是自然堆積下的活性炭試樣經(jīng)過CT掃描后獲得的一幅斷面圖像[58-59]。首先，對(duì)圖像進(jìn)行灰度直方圖均衡化處理，重新分配像素值，使得一定灰度范圍內(nèi)的像素?cái)?shù)量大致相同，增強(qiáng)圖像中孔隙與固體骨架之間的對(duì)比度[57];其次，對(duì)圖像進(jìn)行必要的降噪、濾波和平滑處理，避免圖像生成過程由于噪聲的引入對(duì)辨別孔隙形狀結(jié)構(gòu)造成障礙[60]。經(jīng)過灰度直方圖均衡化處理及均值濾波器平滑處理后，獲得如圖4(b)所示灰度圖;再次，選定閾值，對(duì)圖像進(jìn)行二值化處理，將圖像中每個(gè)像素標(biāo)記為“0”或“1”，其中“0”表示固體，“1”表示孔隙。為了避免圖像邊界處噪聲以及堆積多孔介質(zhì)中壁面附近的邊界效應(yīng)[61]對(duì)分析孔隙結(jié)構(gòu)產(chǎn)生影響，選取圖4(b)中矩形區(qū)域，利用大律算法(Otsu’s method)計(jì)算出背景與目標(biāo)之間方差最大的閾值，對(duì)圖像進(jìn)行閾值分割，獲得如圖4(c)所示的二值圖像。最后，將二值圖片導(dǎo)入MATLAB中得到二維矩陣，二值矩陣中的元素與數(shù)字二值圖像中的像素一一對(duì)應(yīng)，利用cat函數(shù)將一系列二維矩陣重構(gòu)成為三維重構(gòu)矩陣，根據(jù)三維重構(gòu)矩陣構(gòu)建的試樣重構(gòu)模型如圖4(d)所示。

圖4 活性炭試驗(yàn)顯微斷面圖及其圖像處理Fig.4 CT images of an activated carbon sample and corresp- onding processed images

3.2.2表征體元結(jié)構(gòu)提取

生成孔隙結(jié)構(gòu)后，為了計(jì)算試樣孔隙率并提取其三維表征體元，對(duì)三維重構(gòu)矩陣進(jìn)行兩點(diǎn)相關(guān)性分析，多孔介質(zhì)內(nèi)任意兩點(diǎn)的相關(guān)函數(shù)可以表示[62-63]為

S2(r)=

(11)

其中，<>為沿坐標(biāo)軸方向的體積平均;r為多孔介質(zhì)內(nèi)任意兩點(diǎn)的距離;f(x)為一個(gè)結(jié)果為0或1的特征函數(shù)。對(duì)于多孔介質(zhì)來說，當(dāng)x處于孔隙時(shí)，f(x)=1;反之，f(x)=0。兩點(diǎn)相關(guān)函數(shù)可以定量地表征多孔介質(zhì)的孔隙結(jié)構(gòu)，其物理意義可以理解為重構(gòu)矩陣中任意兩個(gè)距離為r的點(diǎn)同時(shí)處于孔隙相的概率。圖5是試樣重構(gòu)矩陣的兩點(diǎn)函數(shù)計(jì)算結(jié)果。當(dāng)r=0時(shí)，兩點(diǎn)相關(guān)函數(shù)計(jì)算結(jié)果表示重構(gòu)矩陣的孔隙率φ(φ=0.45);隨著r逐漸增大，其兩點(diǎn)相關(guān)函數(shù)值快速下降，在r=50時(shí)，兩點(diǎn)相關(guān)函數(shù)值為0.21;當(dāng)r進(jìn)一步增大，相關(guān)函數(shù)計(jì)算結(jié)果基本不變，穩(wěn)定在0.20附近。因此，可以將50像素作為表征體元基本尺寸。當(dāng)基本尺寸大于50時(shí)，表征體元中孔隙結(jié)構(gòu)具有相同的宏觀統(tǒng)計(jì)學(xué)性質(zhì);而當(dāng)尺寸小于50時(shí)，表征體元中孔隙結(jié)果的宏觀特性隨表征體元尺寸大小的改變而變化。

圖5 活性炭試驗(yàn)兩點(diǎn)自相關(guān)函數(shù)計(jì)算結(jié)果Fig.5 Two-point correlation function for the pore phase of the active carbon

4 尺度上聯(lián)方法

從尺度間內(nèi)在關(guān)聯(lián)角度來說，多尺度模擬方法中關(guān)聯(lián)型方法的核心思想是由小尺度描述為上一個(gè)尺度提供本構(gòu)關(guān)系，即不同尺度模型間參量的提煉與傳遞。基于孔隙尺度微觀數(shù)學(xué)模型的計(jì)算結(jié)果，提煉宏觀參數(shù)，獲取滲透率、孔隙率、努謝爾特?cái)?shù)、雷諾數(shù)、普朗特?cái)?shù)等之間的關(guān)系式，可以用于表征體元尺度宏觀數(shù)學(xué)模型的求解。

(1)孔隙率。

煤巖體緩慢升溫過程中孔隙率隨溫度而變化。使用馬弗爐對(duì)特定煤樣進(jìn)行升溫實(shí)驗(yàn)，冷卻后通過電鏡掃描表面方法可獲取局部孔隙率空間分布與溫度的關(guān)系式;或采用工業(yè)CT掃描的數(shù)字重構(gòu)方法，以精確獲取煤巖體升溫后的內(nèi)部結(jié)構(gòu)特征變化，局部孔隙率的變化如式(12)所示。

ε=f(x,y,z,T)

(12)

(2)滲透率和慣性系數(shù)。

對(duì)于特定的孔隙結(jié)構(gòu)，給定不同的壓差條件，先計(jì)算通過煤巖多孔介質(zhì)內(nèi)的流速。再根據(jù)擴(kuò)展達(dá)西定律(式(13))，由多組模擬得到的計(jì)算結(jié)果，利用數(shù)學(xué)方法獲得滲透率K和慣性系數(shù)cF的具體數(shù)值，將這些數(shù)值帶入到式(2)中可以進(jìn)行基于連續(xù)介質(zhì)假定的宏觀滲流場(chǎng)的計(jì)算。如果考慮孔隙率變化對(duì)這兩個(gè)參數(shù)的影響，則通過給定不同的孔隙結(jié)構(gòu)來獲取不同的計(jì)算結(jié)果。煤巖體孔隙結(jié)構(gòu)的構(gòu)造，可以通過人工生成方法或者工業(yè)CT掃描的方法。

(13)

(3)對(duì)流換熱系數(shù)。

針對(duì)特定的孔隙結(jié)構(gòu)，根據(jù)溫度場(chǎng)的模擬計(jì)算結(jié)果，得到流固之間的對(duì)流換熱系數(shù)，用于局部非熱平衡控制方程的計(jì)算。一般來說可通過無量綱準(zhǔn)則數(shù)之間的關(guān)系來構(gòu)建方程，努謝爾特?cái)?shù)可以表示為雷諾數(shù)、普朗特?cái)?shù)和孔隙率的函數(shù)，如式(14)所示:

(14)

5 結(jié) 論

(1)基于煤炭自然發(fā)火的介尺度特性，從數(shù)學(xué)建模角度，提出煤顆粒聚團(tuán)(介尺度Ⅱ)的本質(zhì)就是煤堆的表征體元，明確了顆粒聚團(tuán)的物理邊界，分析了表征體元內(nèi)結(jié)構(gòu)對(duì)于傳熱傳質(zhì)及反應(yīng)過程的影響，實(shí)現(xiàn)了介尺度Ⅱ所在堆積態(tài)煤體層次的研究問題、物理模型和計(jì)算方法在結(jié)構(gòu)上的一致性。

(2)建立了考慮瞬時(shí)孔隙率和高溫輻射換熱特征的表征體元尺度宏觀數(shù)學(xué)模型和考慮內(nèi)部孔隙結(jié)構(gòu)的微觀模型?；诠I(yè)CT掃描和數(shù)值重構(gòu)獲得了煤顆粒聚團(tuán)內(nèi)的真實(shí)孔隙形態(tài)，量化了三維表征體元的尺寸。

(3)考慮表征體元內(nèi)煤巖結(jié)構(gòu)特征變化，利用不同尺度之間參數(shù)的本構(gòu)關(guān)系，提出從孔隙尺度獲取宏觀控制方程中孔隙率、滲透率、慣性系數(shù)和對(duì)流換熱系數(shù)的尺度上聯(lián)方法，從數(shù)學(xué)上實(shí)現(xiàn)了孔隙尺度微觀模型到表征體元尺度宏觀模型之間的信息傳遞，為進(jìn)一步獲得精確的物質(zhì)輸運(yùn)規(guī)律提供前提。