何文浩，郝朝瑜，張亞超，鄧存寶，張 儉

(1.太原理工大學 安全與應急管理工程學院，山西 太原 030024；2.中煤科工西安研究院(集團)有限公司，陜西 西安 710077)

隨著我國礦井進入深部開采階段，瓦斯相關災害(煤與瓦斯突出、瓦斯窒息、瓦斯爆炸等)越來越難以治理和預防，相關文獻統(tǒng)計指出瓦斯相關事故是名副其實的礦井“第一殺手”[1-3]。而其中危害最大的是瓦斯爆炸事故，瓦斯爆炸事故破環(huán)性大和死亡率高，始終是礦井安全管理工作的重點，而抑爆減災技術則是瓦斯爆炸防治工作中的重點內容。

受益于計算機仿真技術的發(fā)展，眾多學者對甲烷氧化鏈式反應進行了細致且精確的分析[4-5]，并結合敏感性分析對反應機理簡化研究[6-8]，大幅降低了瓦斯爆炸抑制理論的難度。自此，引起了各類材料對甲烷爆炸影響的理論研究熱潮。以常見的水對瓦斯爆炸影響為例。最初，陸守香等[9]指出水作為惰性第三體參與自由基的碰撞反應，以物理方式降低自由基濃度，從而抑制瓦斯爆炸。唐建軍[10]通過仿真和試驗結合指出：水不僅有物理抑制作用，還能在一定程度與自由基發(fā)生反應，降低自由基濃度從而抑制爆炸。文獻[11]通過結合量子化學從頭算方法和分子動力學模擬表明：水僅在甲烷氧化鏈式反應初期有抑制作用，而隨著鏈式反應進行到完全激活時，反而會促進甲烷的氧化過程。隨著對反應本質的更深入的剖析，水對瓦斯爆炸的影響才得以漸漸清晰。

硅藻土作為一種綠色環(huán)保，無毒害的抑爆材料[12]，其抑爆效果研究目前主要以試驗方式開展[13-15]，對具體的微觀反應機理則不甚明了。試驗方法來測定抑制效果固然清楚，但難以說明內在抑制機理，不便設計改進方案來增強效果。筆者擬從微觀角度分析硅藻土對瓦斯爆炸的影響，用量化方法表明影響效果，并指出改進方向。

1 計算方法

本文所有反應物、過渡態(tài)和產物結構均在B3LYP-D3/6-311G**水平上進行優(yōu)化，并在同計算級別下開展頻率計算檢驗和內稟反應坐標(IRC)驗證過渡態(tài)正確性。所有反應物和產物均不存在虛頻，過渡態(tài)均只有一個虛頻且聯(lián)通合理反應物和產物?？紤]上述計算級別對能量計算的誤差，本文所有能量計算均在M062X-D3/Def2-TZVP水平進行。此外，在針對HCHO的吸附能計算中，考慮基組的不完備性對弱相互作用體系影響甚大，對其進行基組重疊誤差(BSSE)校正。

1.1 硅藻土模型

硅藻土作為無定形二氧化硅，對其結構的研究較多。文獻[16]建立的無定形二氧化硅理論模型指出，其表面具有3種類型羥基：孤立羥基、連生羥基以及雙生羥基，該模型很好地解釋了硅藻土吸附、催化等各種應用的理論問題。以此理論模型為基礎，無定形二氧化硅的計算化學建模主要分為2種方法：第1種為以DFT理論為代表的小簇模型[17-19]；第2種是以第一性原理為代表的晶面模型[20]。

筆者綜合考慮計算的時長和結果的精度，采用小簇模型來模擬硅藻土表面結構，同時參考關于國產硅藻土結構研究的相關文獻[21]，選擇其中2種常見的硅醇：孤立硅醇(Single Silanol)和連生硅醇(Vicinal Silanol)。模型分子結構如圖1所示，孤立硅醇模型如圖1(a)所示[22]，連生硅醇模型如圖1(b)所示[19](分子結構圖由可視化軟件CYLview導出[23])。

圖1 模型分子結構Fig.1 Molecular structure of models

1.2 硅藻土對瓦斯爆炸影響機理分析

1.2.1 與自由基的反應機理

根據甲烷氧化鏈式反應機理，反應初期自由基濃度對鏈式反應能否傳遞下去有著直接關系[24]，而在眾多自由基中主要以·H，·OH和·HCO為關鍵鏈傳播自由基[7,25-26]。其中，·HCO由穩(wěn)定中間產物HCHO與·H，·OH或其他奪H自由基反應后得到[7,24]。如果能降低·H和·OH數量，則·HCO的數量也會大幅降低。

綜上，筆者著眼于分析硅藻土對以下2個初始鏈式反應的影響。

CH4與·H的反應(亦即表2中的：Ori_H反應)：

CH4與·OH的反應(亦即表2中的：Ori_OH反應)：

若硅藻土與·H和·OH的反應速率能影響甚至超過上述反應，即能有效抑制瓦斯爆炸。由于涉及反應較多，表1為所有可能反應的基本信息。

表1 反應基本信息

表1中反應1～12為硅醇與自由基可能發(fā)生的反應。以反應11,12為例，2者均為連生醇與·OH可能發(fā)生的反應，反應位點F為圖1連生硅醇模型中左側(B為右側硅醇)，產物分別為H2O和H2O2(略去產生的硅醇自由基)，其余反應同理。反應1～4為孤立硅醇與自由基反應，其中反應位點S即代表孤立硅醇。

Reverse反應為反應9產生的連生硅醇自由基與·H的結合反應。

在確定反應過渡態(tài)并驗證其正確性后，進行相關物質的熱力學量計算?？紤]甲烷氧化鏈式反應過程中會引起環(huán)境溫度上升進而影響熱力學量，筆者采用熱力學分析軟件Shermo進行掃描[27]，設定溫度為25～1 000 ℃(即298.15～1 273.15 K)[28]，每25 ℃計算1次，壓力恒定為1 bar(即100 kPa)。由于過渡態(tài)頻率計算結果中存在較多過小頻率(100 cm-1以下)，為提高熱力學量計算精度，計算模式選用Grimme的準RRHO模型[29]?？紤]理論方法的誤差，熱力學頻率校正因子取sclZEP =0.981[30]。

通過熱力學量進一步計算得出反應自由能能壘(計算方法見式(1))，結合過渡態(tài)理論(TST)計算出各反應隨溫度變化的速率常數[31]，比較速率常數的大小來判斷反應相互之間影響大小。

自由能能壘計算公式為

ΔG(T)=GTS(T)-Grea(T)

(1)

式中，GTS(T)為隨溫度變化的過渡態(tài)自由能，kJ/mol；Grea(T)為隨溫度變化的反應物自由能，kJ/mol。

過渡態(tài)理論反應速率常數kTSK計算公式為

(2)

式中,κ為隧道效應校正系數；kB為玻爾茲曼常數，取1.381×10-23J/K；T為絕對溫度，K；h為普朗克常數，取6.626×10-34J·s；R為理想氣體常數，取8.314 47 J/(mol·K)；P0為氣相標準態(tài)壓力；ΔG0為氣相標準態(tài)壓力時，過渡態(tài)的自由能與反應物的自由能之差，kJ/mol。

由于本文反應涉及氫轉移較多，而該類反應有著較強的隧道效應(特別是低溫時)，固對其展開細致計算。

隧道效應校正系數κ計算公式有如下2種情況：

(1)近似的Skodje-Truhlar方法(β<α時可用，對應虛頻小、溫度高時)。

(3)

(4)

(5)

式中，Px為虛頻，cm-1。

(2)完整的Skodje-Truhlar方法(本文用于虛頻超過1 000的反應)。

當β<α時:

(6)

當α<β時:

(7)

式中，Z為正向勢壘，過渡態(tài)減反應物的內能，kJ/mol；V為吸熱反應為產物減反應物內能，kJ/mol,放熱反應取0；NA為阿伏伽德羅常數，取6.02×10-23。

1.2.2 對甲醛的吸附機理

眾多研究表明：在甲烷氧化鏈式反應中，甲醛是不可忽視的關鍵中間體。甲醛的進一步氧化不僅會大幅加速鏈式反應的進程，同時還會釋放大量的熱量來加速反應的進行[4,10,25]。硅藻土由于發(fā)達的微孔結構、豐富的表面極性羥基和低廉的成本，常用作極性有機小分子脫除劑[32-33]，而甲醛正是典型的極性小分子有機物。

反應物分子總是以極性的方式相互接近，進而達到吸附穩(wěn)態(tài)或發(fā)生反應。筆者利用波函數分析軟件Multiwfn[34]，計算分子表面靜電勢，最終通過分子可視化程序VMD呈現(xiàn)并得出可能吸附位點[35]。弱相互作用可視化分析選用IRI方法[36]，同樣通過VMD呈現(xiàn)。

計算弱相互作用時，考慮基組的不完備性，還應加入基組重疊誤差(BSSE)的校正。

考慮基組重疊誤差(BSSE)的吸附能計算公式如下：

E′Binding=EBinding+EBSSE

(8)

EBinding=Ea+b-Ea-Eb

(9)

式中，E′Binding為校正后吸附能，kJ/mol；EBinding為未校正吸附能，kJ/mol；EBSSE為基組重疊誤差校正能，kJ/mol；Ea+b為達到吸附穩(wěn)態(tài)后體系總能量，kJ/mol；Ea，Eb為a，b兩分子單獨穩(wěn)定結構時的能量，kJ/mol。

2 結果與討論

2.1 與自由基的反應機理

從表1可以看出，反應1～12兩兩互為競爭關系，為同種自由基在同一硅醇上的競爭反應。筆者先對這部分反應兩兩之間展開比較分析，判斷優(yōu)勢反應，隨后再與原甲烷氧化鏈式反應比較。

2.1.1 反應過渡態(tài)與能壘

反應1～11過渡態(tài)、虛頻和反應能壘如圖2所示。

以圖2(e)中反應9，10為例，2者均為·OH在連生硅醇B位點上的反應，反應過渡態(tài)結構見上部結構圖。經內稟反應坐標(IRC)驗證，前者反應·OH從B位點抽取氫原子，生成H2O和連生硅醇自由基(缺H)，反應過渡態(tài)虛頻為1 503.638 7 cm-1；后者反應·OH從B位點抽取羥基，生成H2O2和連生硅醇自由基(缺O(jiān)H)，反應過渡態(tài)虛頻為419.717 8 cm-1。反應能壘和能壘差見下部曲線圖，反應10能壘更大(低溫時將近8倍，高溫時10倍左右)，且其隨溫度升高速率明顯大于反應9。其余反應同理。

反應12經反復嘗試搜尋過渡態(tài)和柔性(剛性)掃描均無法發(fā)現(xiàn)合理反應路徑，推測該反應不成立。由于受基團極性影響，·OH在該反應位點附近時，更容易受到B位點H原子吸引，而無法到達F位點合適反應位置。此外，連生硅醇分子間存在氫鍵作用，進一步加大了該反應的難度。綜上，反應11無競爭反應，故圖2(f)中只單獨列出能壘分析。

此外，Ori_H和Ori_OH反應相關文獻中已進行詳細的分析[37]，本文不再贅述，只在后續(xù)分析反應速率常數。

從能壘的角度分析，反應能壘越小，反應物越容易轉變?yōu)檫^渡態(tài)物質，也越容易發(fā)生反應，但其并不是唯一決定反應快慢的量(如隧道效應也會影響反應速率)。嚴格意義上來講，反應速率常數k才是從根本決定反應速率快慢的原因，下面對各反應速率常數進行分析。

2.1.2 反應速率常數

反應速率常數變化如圖3所示。為方便比較各速率之間的影響，圖中以lnk為y軸，以1 000/T為x軸。

在圖3分析中，采用以Δlnk=2.3為影響界限(即范圍在(0.1k,10k))。以Δlnk最小的反應5與反應6為例，2者最接近處位于1 273.15 K(最高溫度)，但劣勢反應仍無法對優(yōu)勢反應產生明顯影響。同理，圖3中的反應1～10中，優(yōu)勢反應為反應1，3，5，7，9。其中反應1，5，7為硅醇與·H反應，反應3，9為硅醇與·OH反應(還應加上未參與比較的反應11)。

下面分別按照與·H或·OH反應的速率，比較上述反應與原鏈式反應Ori_H和Ori_OH速率常數大小。

由圖4可知，2種硅醇與·H的直接反應速率均無法有效影響到原鏈式反應Ori_H。

當環(huán)境溫度隨甲烷氧化放熱升溫到31 ℃時(圖4(b)中1 000/T=3.29 K-1處)，反應9開始對鏈式反應產生影響，即開始減弱鏈式反應作用。隨環(huán)境溫度繼續(xù)升高達到122 ℃時(圖中1 000/T=2.53 K-1處)，反應9的反應速率與Ori_OH相等，即兩者爭奪·OH的速率相等，鏈式反應得以有效減弱。環(huán)境溫度進一步上升到259 ℃時(圖中1 000/T=1.88 K-1處)，反應9速率已為Ori_OH反應的10倍，鏈式反應得到有效抑制，且后續(xù)兩者差距越來越大。此外，在溫度為292 ℃時(圖中1 000/T=1.77 K-1處)，孤立硅醇與·OH反應達到與Ori_OH反應一樣的速率。

圖2 過渡態(tài)與能壘差Fig.2 Transition state and energy barrier difference

圖3 反應速率常數Fig.3 Reaction rate constant

圖4 反應速率常數比較Fig.4 Comparison of reaction rate constants

從前文可知，反應9的產物有2種：H2O和連生硅醇基(缺H)。由于連生硅醇基與·H恰好有合適反應位點，且自由基結合反應的高速特性，猜測Reverse反應成立，并對原鏈式反應有較強影響。通過過渡態(tài)搜尋和剛性掃描，確定該反應無過渡態(tài)。由于反應無過渡態(tài)，速率常數只能利用反應速率公式進行估算(以產物和反應物自由能之差代替能壘，透射系數取1)。估算結果表明：速率最大處為25 ℃時，數量級為10225；速率最小處為1 000 ℃時，數量級為1068；顯然Reaverse反應速率遠大于上述所有反應。另外，由于Reaverse反應依賴于反應9中連生硅醇基的產生速率，即可近似看作兩反應速率相等，反應速率由反應9控制。二者反應相輔相成，反應原理如圖5所示。

圖5 反應流程Fig.5 Reaction flow chart

綜上，硅藻土具有良好抑制甲烷氧化鏈式反應的效果，其中，連生硅醇在低溫時便具備良好的性能?；诖耍裟軐μ烊还柙逋吝M行改性，增加其連生硅醇占比，將會大幅提高抑制效果。

2.2 對甲醛吸附作用

2.2.1 分子表面靜電勢(ESP)分析

從硅醇分子ESP(靜電勢)分布角度出發(fā)，硅醇與甲醛的可能吸附位點如圖6所示。圖中藍色表示靜電勢為負，顏色越深數值越大；紅色表示靜電勢為正，顏色深淺同理。此外，圖中黃色小點表示硅醇分子靜電勢表面極大值點，青色小球表示極小值點。

由圖6可以看出，孤立硅醇模型由于其規(guī)則的結構，僅在羥基側方有與HCHO結合的位點(記為S)。針對連生硅醇模型，以羥基指向左方為正(圖6(c))，分別初步確定了前(F)、后(B)、左(L)、右(R)以及上(U)5個位點。

2.2.2 吸附能計算與分析

各吸附位點構型和相應吸附能如圖7所示。

圖6 硅醇靜電勢極值點分布Fig.6 Distribution on extreme points of silanol electrostatic potential

圖7 吸附構型及能量Fig.7 Adsorption configuration and energy

硅醇分子與HCHO吸附組成了HCHO-硅醇體系，使HCHO的分子運動減弱，降低了HCHO與自由基碰撞的概率，從而抑制后續(xù)鏈式反應的發(fā)生。

通過上表計算表明，連生硅醇的多個吸附位點中，U位點吸附能最大，且明顯大于其他位點。連生硅醇相比孤立硅醇具有更強的吸附性能，更利于脫除鏈式反應中產生的HCHO。

2.2.3 吸附弱相互作用類型分析

利用Multiwfn軟件進行弱相互作用可視化分析，結果如圖8所示(為保證呈現(xiàn)效果，僅導出吸附位點局部圖)。

由色帶對比可知，孤立硅醇與HCHO的吸附作用由3種弱相互作用產生：氫鍵、范德華力和弱氫鍵。氫鍵發(fā)生在HCHO中的O原子與孤立硅醇H原子作用處，此處兩分子間相互吸引最強；弱氫鍵發(fā)生在HCHO中H原子與孤立硅醇Si—O—Si處，吸引力中等；范德華力則夾雜在兩者之間，由于距離各原子較遠，吸引力最小。

連生硅醇與HCHO的吸附只有2種弱相互作用：氫鍵和范德華力。氫鍵發(fā)生在HCHO中的O原子與連生硅醇B位點H原子作用處，此處兩分子間吸引力最大；范德華力發(fā)生在HCHO分子側面相對連生硅醇分子處，吸引力較弱。通過對比原連生硅醇分子內的弱相互作用(圖8(b)，(c))，可知HCHO的加入加強了兩羥基間的氫鍵作用，由發(fā)生吸附前的弱氫鍵作用轉變?yōu)閺姎滏I作用。

圖8 弱相互作用可視化分析Fig.8 Visual analysis of weak interaction

綜上，硅藻土對甲烷氧化鏈式反應中的關鍵中間產物HCHO具有較強的吸附效果。2種硅醇模型的吸附性能差別不大，均與HCHO結合形成氫鍵作用，其中雙生硅醇吸附強度略大一點。

3 結 論

(1)針對甲烷氧化鏈式反應，硅藻土中的連生硅醇結構具備良好的抑制效果。在31 ℃時，便開始影響鏈式反應；在122 ℃時，達到相同速率，有效減弱鏈式反應；在259 ℃及更高溫度時，有效抑制鏈式反應。

(2)硅藻土能與HCHO結合形成氫鍵體系，具備較好的脫除效果，其中連生硅醇在與HCHO的吸附，具有更穩(wěn)定的吸附構型。

(3)不論是活潑的化學性質還是穩(wěn)定的吸附構型，都離不開硅藻土表面的連生硅醇。若設計出經濟高效的改性方法以提高硅藻土中連生硅醇占比，將會大幅增強抑制效果。

(4)參考硅藻土結構，以表面極性羥基更活躍為新型抑制劑的篩選或改性目標。