郝天軒，李 帆，唐一舉,3

(1.河南理工大學(xué)安全科學(xué)與工程學(xué)院，河南 焦作 454000；2.煤炭安全生產(chǎn)河南省協(xié)同創(chuàng)新中心，河南 焦作 454000；3.河南工業(yè)和信息化職業(yè)學(xué)院，河南 焦作 454000)

我國是世界上煤與瓦斯突出最嚴(yán)重的國家之一，具有煤與瓦斯突出危險性的煤礦分布在我國16個省、自治區(qū)和直轄市[1]，因突出事故造成的工人傷亡多次發(fā)生，因此，有效地預(yù)測煤與瓦斯突出將大大提高礦井經(jīng)濟(jì)效益，減少不必要的傷亡。迄今為止，國內(nèi)外學(xué)者已經(jīng)在煤與瓦斯突出機理方面做了許多研究，將其機理主要分為瓦斯主導(dǎo)假說、地壓主導(dǎo)假說、化學(xué)本質(zhì)假說以及綜合作用假說[2]，其中綜合作用假說認(rèn)為煤與瓦斯突出是由地應(yīng)力、煤體中的瓦斯以及煤體自身物理學(xué)性質(zhì)三者綜合作用的結(jié)果，該理論在預(yù)防煤與瓦斯突出的實踐中已經(jīng)得到應(yīng)用。

煤體溫度的變化是地應(yīng)力、瓦斯解吸及解吸后的氣體擴(kuò)散、煤的物理力學(xué)性質(zhì)共同作用的結(jié)果，所以煤體溫度的變化可以作為預(yù)測煤與瓦斯突出的一個有效參數(shù)。國內(nèi)外一些學(xué)者就多物理場耦合、瓦斯吸附解吸和煤體彈性潛能釋放等方面做了相關(guān)研究。趙陽升等[3]和趙延林等[4]介紹了在鹽礦開采、高溫巖體地?zé)衢_采等方面利用固流熱耦合的模型進(jìn)行數(shù)值模擬；陶云奇[5-6]構(gòu)建了含瓦斯煤的THM耦合模型，該模型實現(xiàn)了含瓦斯煤的雙向完全耦合；牛國慶等[7]通過煤體吸附、解吸CO2和N2的實驗驗證了煤體吸附瓦斯過程是放熱過程，瓦斯解吸過程是吸熱過程，溫度變化幅度隨壓力變化幅度增加而增加；郭立穩(wěn)等[8-9]通過理論分析結(jié)合實驗認(rèn)為，煤體溫度在突出過程中升高是由地應(yīng)力破碎煤體，使得彈性勢能釋放造成的溫度升高，而溫度降低是由于瓦斯解吸擴(kuò)散引起的，變化是先升高后降低的連續(xù)變化；前人有通過多物理場的耦合來模擬煤巖或其他礦類開采，也有通過實驗或理論介紹煤體溫度變化的影響因素，其中都有涉及溫度變化，但影響階段和環(huán)境不盡相同。本文重點在于結(jié)合影響煤體溫度多種因素共同作用的結(jié)果，研究含瓦斯煤體在不同瓦斯壓力下單軸受壓過程中表面溫度變化規(guī)律。

影響煤與瓦斯突出的應(yīng)力條件十分復(fù)雜，本文通過COMSOL軟件結(jié)合相關(guān)學(xué)科知識構(gòu)建熱流固耦合模型，模擬煤體在單軸受壓過程中的溫度變化并研究溫度變化規(guī)律，能夠較客觀地描述溫度場變化，由此可見煤體溫度參數(shù)可以用來反映煤體受壓情況，在利用煤體溫度預(yù)測預(yù)防煤與瓦斯突出的研究上提供一定理論依據(jù)，為預(yù)測預(yù)防煤與瓦斯突出工作提供幫助。

1 模型建立

1.1 基本假設(shè)

由于影響成煤的地質(zhì)因素較多且復(fù)雜，使得煤體呈現(xiàn)非均質(zhì)性，煤體中的瓦斯?jié)B流運移也難以控制，同時在含瓦斯煤的孕災(zāi)過程中，沒有直接熱源作用，所以對溫度造成影響的主要因素即應(yīng)力場作用、滲流作用、瓦斯解吸吸熱效應(yīng)以及煤體內(nèi)氣體擴(kuò)散的吸熱效應(yīng)等。 考慮到以上條件作出如下假設(shè)：①煤與瓦斯突出由地應(yīng)力、瓦斯和煤體自身物理學(xué)性質(zhì)綜合作用，其中地應(yīng)力起主導(dǎo)作用；②含瓦斯煤體為均質(zhì)且各向同性的線彈性體；③含瓦斯煤骨架的有效應(yīng)力變化遵循Terzaghi有效應(yīng)力規(guī)律；④瓦斯在煤層中的滲流規(guī)律符合Darcy’s Law[10]；⑤將瓦斯視為理想氣體，用氣體狀態(tài)方程表示其密度；⑥煤體中吸附狀態(tài)瓦斯服從Langmuir吸附平衡方程；⑦煤體的變形是微小的，煤體處于線彈性變形階段，遵循廣義胡克定律；⑧溫度對孔隙率以及應(yīng)變的影響較小，因此暫不考慮溫度對孔隙率和應(yīng)變的影響。

1.2 流-固-熱耦合模型

流-固-熱耦合模型由滲流場、應(yīng)力場以及溫度場三個物理場的控制方程組成，分別控制滲流場、應(yīng)力場和溫度場的變化。

1) 含瓦斯煤應(yīng)力場控制方程。依據(jù)假設(shè)，認(rèn)為煤體受壓過程為彈性變形，符合胡克定律，結(jié)合有效應(yīng)力理論，含瓦斯煤應(yīng)力場控制方程見式(1)[11]。

(i，j=1,2,3)

(1)

式中：G為剪切彈性模量；u為位移；Fi為體積力；α為Boit系數(shù)；P為氣體壓力。

由于溫度對應(yīng)變影響較小，所以該式暫不考慮溫度對應(yīng)變的影響，主要體現(xiàn)瓦斯壓力與地應(yīng)力對煤體造成的變形控制。

2) 含瓦斯煤滲流場控制方程。根據(jù)假設(shè)，瓦斯在煤體中滲流符合Darcy’s Law，吸附狀態(tài)瓦斯符合Langmuir吸附方程，再根據(jù)質(zhì)量守恒以及公式變形得到滲流場控制方程式，見式(2)[12]。

(2)

式中：φ為孔隙率；Ks為煤體體積模量；k為滲透率；μ為黏滯系數(shù)；P為煤體中瓦斯壓力梯度；T為氣體溫度；Qd為含瓦斯煤的瓦斯吸附量，由于不考慮溫度對孔隙率的影響，φ=φ(P,e)。

式(2)中左側(cè)第一項為含瓦斯煤中瓦斯壓力的影響項；第二項為溫度對滲流場的影響項；第三項為應(yīng)變對滲流場的影響項；第四項為含瓦斯煤所吸附瓦斯量對滲流場的影響。

3) 含瓦斯煤溫度場控制方程。依據(jù)熱力學(xué)第一定律，系統(tǒng)內(nèi)能的增加等于系統(tǒng)吸收熱量和對系統(tǒng)做功的總和，其中系統(tǒng)吸熱Q為正值，系統(tǒng)放熱Q為負(fù)值，同理對系統(tǒng)做功W為正值，系統(tǒng)對外做功W為負(fù)值，由此得到溫度場控制方程，見式(3)。

(3)

式中：φ為孔隙率；ρs為煤體密度；cs為煤體的比熱容；ρl為瓦斯密度；cl為瓦斯的比熱容；kt為熱傳導(dǎo)系數(shù)；其中Q=Q(P)(為負(fù))，W=W(σ，u)(為正)，左側(cè)第一項為含瓦斯煤系統(tǒng)內(nèi)的內(nèi)能變化，第二項為瓦斯壓力與應(yīng)變對溫度場的影響項，第三項為熱傳導(dǎo)項。

聯(lián)立式(1)、式(2)和式(3)，即為含瓦斯煤的熱流固耦合模型，見式(4)。

(4)

由于孔隙率和滲透率也受瓦斯壓力和滲透率的影響，所以除了上述的耦合模型，還需聯(lián)立孔隙率和滲透率模型才可以進(jìn)行模擬計算，孔隙率和滲透率模型見式(5)。

(5)

1.3 幾何模型建立

根據(jù)上述的假設(shè)條件以及含瓦斯煤體的相關(guān)參數(shù)，通過模擬煤樣進(jìn)行單軸加壓試驗來研究含瓦斯煤體在外力受壓過程中的溫度變化規(guī)律，建立含瓦斯煤體的幾何模型(圖1)。 模型為圓柱形煤柱，構(gòu)建時選擇二維旋轉(zhuǎn)體。 圖1為矩形旋轉(zhuǎn)面，圖2為旋轉(zhuǎn)得到的煤體模型，旋轉(zhuǎn)面寬25 mm，模型高100 mm。在COMSOL軟件中選擇多孔介質(zhì)傳熱、固體力學(xué)以及達(dá)西定律三個模塊，其中固體力學(xué)模塊中選用線彈性模塊，縱向破壞變形比例設(shè)置為5%，變形速率按照彈性變形進(jìn)行公式編輯。根據(jù)單軸壓裂煤體試驗，設(shè)置模擬時間為15 min，模型底部為固定邊界，側(cè)表面為自由邊界，頂部設(shè)定指定的位移和荷載，整體上設(shè)定煤體載荷。由于環(huán)境溫度的影響不可控并且對模擬試驗影響較大，所以該模擬試驗不考慮含瓦斯煤體系統(tǒng)之外的熱源影響，設(shè)定煤體初始溫度為293.15 K(20 ℃)，同時設(shè)定含瓦斯煤體的初始瓦斯壓力為1 MPa，自由邊界上壓力為1 atm，模型主要參數(shù)見表1。

圖1 矩形旋轉(zhuǎn)面Fig.1 Rectangular rotating face

圖2 含瓦斯煤體煤柱立體圖Fig.2 Stereogram of coal pillar with gas

表1 含瓦斯煤體的幾何模型相關(guān)參數(shù)Table 1 Related parameters of geometric model ofcoal body with gas

1.4 定解條件

在解決數(shù)學(xué)物理問題時，通常需要(偏)微分方程對問題進(jìn)行描述，而(偏)微分方程一般都有無窮多個解，但具體的物理問題只有唯一解或特定解。該模擬實驗的數(shù)學(xué)模型為偏微分方程組，包括應(yīng)力場方程、滲流場方程以及溫度方程，三個方程相互耦合，想要求解方程組，需要給出定解條件即初始條件和邊界條件。

邊界條件分為第一類邊界條件、第二類邊界條件和第三類邊界條件。第一類邊界條件用以直接描述物理系統(tǒng)邊界上待求解的物理量；第二類邊界條件描述物理系統(tǒng)邊界上物理量的倒數(shù)情況；第三類邊界條件可看做第一類和第二類的線性加和。如下所述給出模型的定解條件。

1) 應(yīng)力場的定解條件。 應(yīng)力場邊界條件見式(5)。

(5)

應(yīng)力場初始條件見式(6)。

u|t=0=u0

(6)

2) 滲流場的定解條件。滲流場的邊界條件見式(7)。

Ps=const

(7)

滲流場初始條件見式(8)。

P|t=0=P0

(8)

3) 溫度場的定解條件。 溫度場初始條件見式(9)。

T|t=0=T0

(9)

2 模擬結(jié)果及分析

如圖3和圖4所示，含瓦斯煤體在加壓1.0 min后，整體溫度開始降低，降低幅度在0.100 ℃左右，煤體表面與煤體中心的溫度變化差異不大；含瓦斯煤體在加壓2.5 min后，受壓煤體中心升溫趨勢影響，使得煤體中心內(nèi)部溫度下降速度、幅度都有一定程度的減緩，但煤體下半部表面仍由于瓦斯的解吸、擴(kuò)散作用溫度加速下降；含瓦斯煤體在加壓3.1 min后，在煤體內(nèi)部的底部形成一個較低溫區(qū)域，該區(qū)域與煤體溫度較高的區(qū)域溫度差越來越大，剛出現(xiàn)該區(qū)域時，溫差僅有0.002 ℃，隨著煤體的受壓變形，溫差達(dá)到了0.030 ℃，在模擬試驗的最后，煤體整體溫度差由-0.113 ℃降低到了-1.560 ℃。

圖3 1.0～3.1 min溫度分布云圖Fig.3 Temperature distribution cloud chart from 1.0 min to 3.1 min

圖4 9.0～15.0 min溫度分布云圖Fig.4 Temperature distribution cloud chart from 9.0 min to 15.0 min

圖5～圖7是選取模擬含瓦斯煤體煤柱側(cè)表面中不同坐標(biāo)點的溫度變化曲線圖，圖8為三個坐標(biāo)點溫度曲線的集合，可以看出含瓦斯煤體溫度隨著煤體受壓持續(xù)下降。由圖5～圖7可知系統(tǒng)溫度下降幅度最大為1.55 ℃，在煤體表面溫度下降最低的坐標(biāo)為點(25，15)，下降幅度為1.489 ℃，點(25，0)和點(25，55)溫度下降分別為1.472 ℃和1.474 ℃。在此過程中，0～5.2 min階段，溫度下降了0.60 ℃；5.2～10.1 min階段，溫度下降了0.50 ℃；10.1～15.0 min階段，溫度下降了0.45 ℃。如圖8所示，在6.0 min左右時，點(25，55)的溫度曲線與點(25，15)的溫度曲線較為接近；在10.6 min左右時，三個點的溫度呈近乎平行的趨勢；在12.0 min左右時，點(25，55)與點(25，0)的溫度曲線開始靠近，直到最后結(jié)束。含瓦斯煤體表面溫度下降的幅度隨受壓時間延長而逐漸降低，不同位置的溫度差異較小，煤體的溫度整體呈下降趨勢。

圖5 含瓦斯煤體側(cè)表面坐標(biāo)為(25，0)處的溫度變化曲線圖Fig.5 Temperature change curve with (25,0)coordinate of gas-containing coal side surface

圖6 含瓦斯煤體側(cè)表面坐標(biāo)為(25，15)處的溫度變化曲線圖Fig.6 Temperature change curve with (25,15)coordinate of gas-containing coal side surface

圖7 含瓦斯煤體側(cè)表面坐標(biāo)為(25，55)處的溫度變化曲線圖Fig.7 Temperature change curve with (25,55)coordinate of gas-containing coal side surface

圖8 含瓦斯煤體側(cè)表面溫度變化曲線圖Fig.8 Surface temperature change curve ofgas-containing coal side

圖9為初始瓦斯壓力分別為0.1 MPa、0.5 MPa、1.0 MPa、1.5 MPa時選取煤體矩形旋轉(zhuǎn)面坐標(biāo)為點(25，15)處的溫度變化曲線。當(dāng)煤體內(nèi)瓦斯壓力為0.1 MPa時，煤體受壓后表面溫度呈上升趨勢，上升幅度較小為0.2 ℃左右。由此可見，含瓦斯煤體中瓦斯初始壓力越大即瓦斯含量越多，在單軸受壓時溫度下降幅度也隨之增加，其中當(dāng)煤體的初始瓦斯壓力為1.5 MPa時，溫度下降幅度在2.4 ℃左右，當(dāng)煤體中含有瓦斯時，瓦斯壓力越小，溫度下降幅度隨時間變化而降低。

圖9 初始瓦斯壓力不同的煤體表面溫度曲線Fig.9 Surface temperature curve of coal withdifferent initial gas pressure

結(jié)合以上的模擬數(shù)據(jù)進(jìn)行分析，含瓦斯煤體受壓初期，煤體表面瓦斯快速解吸，引起瓦斯擴(kuò)散吸收熱量，且瓦斯解吸量和氣體擴(kuò)散量較大，吸收能量較多，與外界能量交換速率快，所以降溫幅度較大，使得煤體在初期時表面溫度比煤體內(nèi)部溫度較低。一方面隨著外力作用，煤體被壓縮，外力做功釋放能量，在一定程度上減緩了溫度降低速度；另一方面瓦斯的解吸、擴(kuò)散隨時間推移使得煤體內(nèi)吸附瓦斯減少，所以后期降溫過程沒有模擬開始時的速率快，降溫幅度也有所減小。在受壓過程中，煤體下半部分為應(yīng)力集中區(qū)，前期解吸的瓦斯主要從這些部位滲流出煤體，所以溫度降低較為明顯，同時在這些部位容易出現(xiàn)煤體摩擦擠壓，導(dǎo)致煤體局部溫度有上升趨勢。

3 結(jié) 論

1) 瓦斯解吸和瓦斯擴(kuò)散吸收熱量導(dǎo)致煤體溫度降低，含瓦斯煤體系統(tǒng)外部受力是系統(tǒng)溫度升高的主要原因，導(dǎo)致煤體的彈性潛能釋放溫度升高，其中溫度變化幅度在0.2 ℃左右。瓦斯解吸、擴(kuò)散吸收的熱量大于煤體受外力作用釋放的熱量。

2) 根據(jù)多個物理場耦合結(jié)果，煤體在受壓過程中，含瓦斯煤體溫度整體降低，其中底部降低幅度更為明顯。煤體中瓦斯壓力為0.5 MPa、1.0 MPa、1.5 MPa時降溫幅度分別為0.58 ℃、1.47 ℃、2.40 ℃，隨著瓦斯壓力增加，溫度的變化幅度逐漸升高，且溫度降低幅度隨著煤體受壓時間越來越小，當(dāng)瓦斯壓力不大時表現(xiàn)的更為明顯。

3) 本文主要研究煤體中所含瓦斯壓力的不同對煤體溫度的影響，今后應(yīng)繼續(xù)研究在此過程中對煤體的其他因素如孔隙率、煤的堅固性系數(shù)等因素的影響，完善相關(guān)理論的研究。