劉殿德 陳志峰 糜行軍

（滕州郭莊礦業(yè)有限責(zé)任公司，山東 滕州 277500）

探明含水層間水力聯(lián)系是地下水資源和水文地質(zhì)勘探研究的基礎(chǔ)[1]，是預(yù)防礦井水害的重要工作之一[2]。含水層之間水力聯(lián)系是一個(gè)復(fù)雜的系統(tǒng)工程[3]，需要綜合多因素多方法確定，包括含水層水化學(xué)分析、含水層水位標(biāo)高變化分析、構(gòu)造賦存情況等因素。對(duì)于水化學(xué)分析方法，對(duì)比不同含水層的地下水礦物離子成分及其變化，有時(shí)會(huì)結(jié)合現(xiàn)場(chǎng)示蹤試驗(yàn)方法，來判斷含水層之間的水力聯(lián)系[4]；對(duì)于含水層水位標(biāo)高變化分析方法，可以對(duì)含水層進(jìn)行抽水試驗(yàn)，研究不同含水層之間的水頭高度變化，判斷同一含水層不同區(qū)域以及不同含水層之間的水力聯(lián)系。

在一些工程現(xiàn)場(chǎng)進(jìn)行抽水試驗(yàn)時(shí)，通常還會(huì)利用水溫變化判斷含水層之間的水力聯(lián)系。由于地溫梯度的影響，含水層之間的溫度也會(huì)不同，通常情況下深度越大，溫度越高，可以利用溫度來判斷地下水所處的層位。但是，在含水層水力聯(lián)系研究方面，從理論上對(duì)地溫場(chǎng)的變化規(guī)律研究較少，缺少相關(guān)物理力學(xué)模型。針對(duì)上述存在的問題，基于多孔介質(zhì)滲流、熱力學(xué)相關(guān)理論，建立溫度-滲流耦合條件下含水層水力聯(lián)系數(shù)值模擬模型，研究不同含水層水力聯(lián)系狀態(tài)、不同回采距離條件下滲流場(chǎng)與溫度場(chǎng)的變化規(guī)律。

1 控制方程

1856 年，法國水力學(xué)家達(dá)西（H.Darcy）通過大量的實(shí)驗(yàn)，得到線性滲透定律，即著名的達(dá)西定律：

式中：Q為通過滲流流量，即實(shí)驗(yàn)過程中通過砂柱各斷面的流量，m3/s；K為滲透系數(shù)，m/s；A為過水?dāng)嗝?，m2；I為水力梯度，等于過水?dāng)嗝嬷g的距離與水頭損失的比值?？梢钥闯?，滲透系數(shù)可以理解為水力梯度為1 時(shí)（單位水力梯度）的滲流速度。同時(shí)應(yīng)注意到達(dá)西定律實(shí)驗(yàn)得出的是滲透系數(shù)，而不是滲透率。

上述達(dá)西定律的形式還不能夠進(jìn)行數(shù)值模擬計(jì)算，需要將其轉(zhuǎn)換成微分形式。達(dá)西定律也可以用壓力P和流體的密度ρf表示，即可以寫成如下形式[5]：

式中：K′即為滲透性常數(shù)，m3·s/kg；ρf為地下水的密度，kg/m3；P為水的壓力，Pa；h為滲流斷面之間的高度差，m。可以看出，其量綱包含了質(zhì)量、長度和時(shí)間單位，表示了某種介質(zhì)對(duì)某特定的流體的滲透能力，它的大小由介質(zhì)和流體兩者的性質(zhì)而定[5]。但是，這種形式還是無法將多孔介質(zhì)的影響與流體的影響分開，于是，便進(jìn)一步引入一個(gè)新的物理量，其除以流體的黏度的值等于滲透性常數(shù)，即：

式中：κ為多孔介質(zhì)的滲透率，m2；μ為地下水的動(dòng)力黏度，Pa·s。滲透率的性質(zhì)取決于多孔介質(zhì)的內(nèi)部結(jié)構(gòu)，與流體的性質(zhì)無關(guān)。

接下來將式（3）中的滲透性常數(shù)替換為滲透率，同時(shí)將壓力差與h的比值表示為壓力梯度，就可以將達(dá)西定律表示為常見的微分形式，即：

式中：v為地下水的滲流速度，m/s；p為孔隙水壓力，Pa；g為重力加速度，m/s2。

地下水在地下巖體中的流動(dòng)還需要滿足連續(xù)性方程，即：

上式用來描述地下水在含水層中的流動(dòng)問題，多孔介質(zhì)內(nèi)傳熱問題通常用基于局部熱平衡理論和體積平均方法導(dǎo)出的多孔介質(zhì)傳熱方程描述，即[6]：

式中：cf為地下水的恒壓比熱容，J/（kg·K）；cs為固體的恒壓比熱容，J/（kg·K）；qeff為多孔介質(zhì)導(dǎo)熱通量，J/（s·m2）；（ρcp）eff為等效單位體積熱容，J/（m3·K）；ks為固體的導(dǎo)熱系數(shù)，W/（m·K）；kf為流體的導(dǎo)熱系數(shù)，W/（m·K）；keff為多孔介質(zhì)的等效導(dǎo)熱系數(shù)，W/（m·K）；T為溫度，K。對(duì)于多孔介質(zhì)，固體骨架和流體同時(shí)存在于一個(gè)體積空間內(nèi)，兩者熱力學(xué)特征存在差異，需要使用等效的熱力學(xué)參數(shù)將兩者統(tǒng)一起來，得到多孔介質(zhì)的能量方程。在累積項(xiàng)中，多孔介質(zhì)內(nèi)的等效單位體積熱容可以使用式（10）進(jìn)行計(jì)算。

2 數(shù)值模型建立

多物理場(chǎng)耦合數(shù)值模擬軟件COMSOL Multiphysics 能夠?qū)⒉煌锢韴?chǎng)的控制方程耦合，模擬研究多物理場(chǎng)耦合問題。使用COMSOL Multiphysics 多孔介質(zhì)流動(dòng)模塊和多孔介質(zhì)傳熱模塊，將達(dá)西定律、多孔介質(zhì)流動(dòng)控制方程以及多孔介質(zhì)傳熱方程耦合，建立含水層水力聯(lián)系溫度場(chǎng)-滲流場(chǎng)耦合模型，實(shí)現(xiàn)地下水在地下巖體多孔介質(zhì)中的流動(dòng)過程和傳熱過程的耦合。具體包括建立幾何模型、設(shè)置物理場(chǎng)與控制方程、設(shè)置初始與邊界條件以及制定研究方案求解溫度場(chǎng)-滲流場(chǎng)耦合控制方程。

軟件中的達(dá)西定律接口用于描述流體在完全飽和的多孔介質(zhì)中通過間隙的流動(dòng)，這種運(yùn)動(dòng)由壓力梯度驅(qū)動(dòng)，流體的剪切應(yīng)力引起的動(dòng)量傳遞可以忽略不計(jì)?？梢允褂眠_(dá)西定律接口計(jì)算壓力，然后根據(jù)壓力梯度、流體黏度和滲透率來確定速度場(chǎng)。其控制方程就是式（6）。

軟件中的多孔介質(zhì)傳熱接口是通過傳導(dǎo)、對(duì)流和彌散進(jìn)行的。彌散是由液體在多孔介質(zhì)中的曲折路徑引起的，并通過包含平均對(duì)流以外的其他效應(yīng)進(jìn)行描述。在許多情況下，固相可以由電導(dǎo)率不同的多種材料組成，也可能會(huì)涉及許多不同的流體。多孔介質(zhì)傳熱接口可以自動(dòng)分析這些因素，并提供用于計(jì)算有效傳熱屬性的混合規(guī)則?？紫犊臻g中的流體也可以經(jīng)歷一個(gè)或多個(gè)相變，這在模擬凍土過程時(shí)很有意義?？梢允褂脤ｉT的相變材料特征通過指定兩種材料和多個(gè)相變屬性（相變溫度、轉(zhuǎn)變間隔和潛熱）來模擬此過程和類似的過程。其控制方程就是式（8），其參數(shù)可以在相應(yīng)的節(jié)點(diǎn)處定義。

2.1 幾何模型

COMSOL Multiphysics 軟件包提供了豐富的幾何建模工具，支持使用實(shí)體對(duì)象、表面、曲線和布爾操作來創(chuàng)建零件?？梢酝ㄟ^操作序列來定義幾何，序列中的每個(gè)操作都可以接收輸入?yún)?shù)，方便在多物理場(chǎng)模型中輕松進(jìn)行編輯和參數(shù)化研究。幾何與定義的物理場(chǎng)設(shè)置之間是完全相互關(guān)聯(lián)的，只要幾何發(fā)生變化，系統(tǒng)便會(huì)自動(dòng)將相關(guān)變化傳遞到整個(gè)關(guān)聯(lián)的模型設(shè)置中?？梢詫⒉牧嫌蚝捅砻娴葞缀螌?shí)體進(jìn)行分組，創(chuàng)建不同的選擇，并在定義物理場(chǎng)、網(wǎng)格劃分以及繪圖等后續(xù)操作中使用這些選擇。不僅如此，還可以使用一系列操作來創(chuàng)建參數(shù)化幾何零件（包括相關(guān)選擇），然后將它們存儲(chǔ)到“零件庫”中，以便在多個(gè)模型中重復(fù)使用。

為了研究不同回采距離下含水層之間的水力聯(lián)系以及滲流場(chǎng)、溫度場(chǎng)之間的關(guān)系，建立三維幾何模型。模型長500 m、高400 m，沿煤層走向厚度600 m，煤層厚度5 m，工作面寬150 m（見圖1）。

圖1 三維幾何模型

2.2 物理場(chǎng)設(shè)置

COMSOL 軟件提供了一系列預(yù)定義的物理場(chǎng)接口，用于模擬各種物理現(xiàn)象，包括許多常見的由多個(gè)物理場(chǎng)共同作用引起的現(xiàn)象。每個(gè)物理場(chǎng)接口都提供專用于相關(guān)科學(xué)或工程領(lǐng)域的特定設(shè)置。當(dāng)選擇所需的物理場(chǎng)接口之后，軟件會(huì)給出適用的研究類型建議，例如瞬態(tài)或穩(wěn)態(tài)求解器。選定求解器后，將對(duì)數(shù)學(xué)模型、求解器序列以及可視化和結(jié)果設(shè)置進(jìn)行適當(dāng)?shù)臄?shù)值離散化。同時(shí)，可以對(duì)所有設(shè)置進(jìn)行編輯操作。

滲流場(chǎng)選擇達(dá)西定律，包含了達(dá)西定律方程和連續(xù)性方程，在其內(nèi)部各節(jié)點(diǎn)可以設(shè)置流體的密度、動(dòng)力黏度等參數(shù)，同時(shí)可以設(shè)置不同域?qū)?yīng)的控制方程，本模型所有域均為達(dá)西定律方程；溫度場(chǎng)選擇多孔介質(zhì)傳熱方程，在其內(nèi)部各節(jié)點(diǎn)可以設(shè)置流體和固體的熱力學(xué)參數(shù)，比如導(dǎo)熱系數(shù)、比熱容，并選擇有效導(dǎo)熱系數(shù)計(jì)算多孔介質(zhì)的熱力學(xué)參數(shù)，同時(shí)本模型所有域均設(shè)置為多孔介質(zhì)傳熱控制方程。

2.3 初始條件與邊界條件

分別對(duì)滲流場(chǎng)和溫度場(chǎng)設(shè)置初始條件和邊界條件。滲流場(chǎng)的邊界條件設(shè)置分為兩種情況：一種情況是設(shè)置垂直方向的壓力梯度，模型上部邊界為定壓力邊界，水壓值根據(jù)研究方案確定，工作面邊界設(shè)置為0.1 MPa，其余邊界設(shè)置為無流動(dòng)；另一種情況是設(shè)置水平方向的壓力梯度，給定模型左側(cè)邊界為定壓力邊界，水壓值根據(jù)研究方案確定，模型右側(cè)邊界設(shè)置為0.1 MPa，其余邊界設(shè)置為無流動(dòng)。對(duì)于溫度場(chǎng)，整個(gè)區(qū)域的地溫梯度設(shè)置為0.025 ℃/m，模型的上下邊界設(shè)置為定溫度邊界，分別為25 ℃和35 ℃。

2.4 研究方案

COMSOL Multiphysics 軟件可以提供瞬態(tài)、穩(wěn)態(tài)或特征頻率研究，同時(shí)可以進(jìn)行掃描、優(yōu)化和估計(jì)研究，通過參數(shù)化掃描來運(yùn)行，其中可以包含模型中的一個(gè)或多個(gè)參數(shù)，涉及幾何參數(shù)、物理場(chǎng)定義中的設(shè)置等。可以使用不同的材料及其定義的屬性來執(zhí)行掃描，也可以對(duì)一組定義的函數(shù)執(zhí)行掃描。

本模型采用穩(wěn)態(tài)研究，為了研究不同回采距離范圍內(nèi)、不同水壓、不同壓力梯度方向條件下滲流場(chǎng)和溫度場(chǎng)的空間分布，對(duì)開采距離、水壓值進(jìn)行參數(shù)化掃描，兩種情況下滲流場(chǎng)上游水壓設(shè)置為0.5～3 MPa。

3 結(jié)果

不同工作面回采距離會(huì)影響工作面圍巖溫度場(chǎng)的分布，不同水壓值和壓力梯度方向決定了滲流速度的變化，地下水滲透壓和滲流速度會(huì)改變地下巖體多孔介質(zhì)的熱力學(xué)性質(zhì)，進(jìn)而決定了地溫的變化。

3.1 垂直方向溫度場(chǎng)分布

圖2 展示了回采350 m 時(shí)，沿回采方向不同位置處、不同水壓下溫度場(chǎng)的分布，沿回采方向（Z軸方向）選取3 個(gè)截面，分別為Z=0 m、Z=300 m和Z=600 m。Z=600 m 截面表示未回采區(qū)域，可以看出，水平方向溫度相同；Z=0 m 和Z=300 m 表示已回采區(qū)域的不同位置，Z=0 m 表示切眼位置，可以看出，工作面上方一定范圍內(nèi)溫度低于兩側(cè)溫度，相同區(qū)域內(nèi)Z=0 m 截面的溫度低于Z=300 m 截面的溫度。

圖2 垂直方向壓力梯度條件下溫度場(chǎng)分布

為了更清晰地展示沿回采方向不同位置處工作面上方的溫度變化規(guī)律，選取工作面上方附近與Z軸平行的一條測(cè)線，其端點(diǎn)坐標(biāo)為（300, 160, 0）和（300, 160, 600），測(cè)線上的溫度變化見圖3?？梢钥闯觯畨河?.5 MPa 增加至3 MPa 過程中，在切眼處（Z=0 m）溫度從29.2 ℃降至25.8 ℃，在Z=600 m 處溫度從29.4 ℃降至26.5 ℃，在Z 軸方向上距離切眼越遠(yuǎn)，工作面上方相同區(qū)域內(nèi)溫度越高，同時(shí)壓力越大，工作面附近溫度越低。產(chǎn)生以上現(xiàn)象的原因是，壓力梯度的方向是垂直方向，地下水在含水層之間的聯(lián)系也是主要發(fā)生在垂直方向，工作面上部區(qū)域的地下水補(bǔ)給主要來自上部含水層，上部含水層由于埋深較淺，地溫較低，地下水由含水層流入工作面的過程中會(huì)與周圍巖體發(fā)生熱量交換，降低工作面附近的溫度。由上述現(xiàn)象可以得出，若工作面附近溫度降低，則反映了地下水來自賦存較淺的含水層，上覆淺部含水層和深部含水層之間發(fā)生水力聯(lián)系。同時(shí)，對(duì)工作面上部區(qū)域進(jìn)行溫度監(jiān)測(cè)，若工作面上部同一高度位置溫度發(fā)生明顯的溫度變化，則說明該區(qū)域的滲流場(chǎng)發(fā)生變化，需要加強(qiáng)對(duì)水害的監(jiān)測(cè)預(yù)警。

圖3 不同水壓下測(cè)線上的溫度分布

3.2 水平方向滲流速度及溫度場(chǎng)分布

圖4 展示了回采350 m 時(shí)，沿回采方向不同位置處、不同水壓下溫度場(chǎng)的分布。沿回采方向同樣選取上述3 個(gè)截面，分別為Z=0 m、Z=300 m 和Z=600 m?？梢钥闯?，3 個(gè)截面表示的區(qū)域溫度分布相同，沒有變化。同時(shí)，在不同水壓和回采方向上，測(cè)線上的溫度也沒有發(fā)生變化。產(chǎn)生這種現(xiàn)象的原因是，壓力梯度的方向是水平方向，地下水在含水層之間沒有水力聯(lián)系，水力聯(lián)系主要發(fā)生在同一巖層內(nèi)，地下水補(bǔ)給主要來自同一含水層，同一含水層在水平方向上由于所處的深度相同，溫度也相同，地下水與其在水平方向上流動(dòng)所經(jīng)過的區(qū)域溫度相同，不會(huì)發(fā)生熱量交換，所以不會(huì)使流經(jīng)區(qū)域的溫度升高或者降低，因而，溫度場(chǎng)不會(huì)發(fā)生變化。由上述現(xiàn)象可以得出，含水層之間沒有水力聯(lián)系，水力聯(lián)系只發(fā)生在水平方向，不會(huì)引起溫度場(chǎng)的變化。

圖4 水平方向壓力梯度條件下溫度場(chǎng)分布

4 結(jié)論

1）溫度場(chǎng)-滲流場(chǎng)耦合數(shù)值模擬模型能夠模擬回采方向上不同位置、不同含水層水力聯(lián)系情況下的溫度場(chǎng)分布，溫度場(chǎng)的變化規(guī)律能夠用來判斷含水層水力聯(lián)系規(guī)律。

2）當(dāng)垂直方向上含水層之間發(fā)生水力聯(lián)系時(shí)，工作面附近溫度會(huì)降低，且距離開切眼越遠(yuǎn)，溫度越高；當(dāng)?shù)叵滤诤畬又g沒有水力聯(lián)系時(shí)，在同一含水層不同區(qū)域之間有水平方向的水力聯(lián)系時(shí)，不會(huì)引起工作面附近溫度場(chǎng)的變化。

3）工作面附近溫度場(chǎng)的變化能夠反映含水層之間的水力聯(lián)系情況，可以用來判斷頂板含水層不同區(qū)域之間以及不同含水層之間是否發(fā)生水力聯(lián)系，進(jìn)而可以用來預(yù)防礦井水害。