基于COMSOL Multi-physic的瓦斯抽采地面井的流場數值分析

翟清偉

（兗州煤業股份有限公司鮑店煤礦，山東 濟寧 273500）

0 引 言

COMSOL Multi-physic是一種多重物理量耦合軟件，該軟件是由MATLAB軟件工具箱發展而來的，以有限元方法進行分析求解，其優點在于高度的靈活性,強大的求解器和較高的計算精度，進行求解時只需要將所建立的數學模型輸入軟件的PDE模塊中，設定求解域，指定邊界條件并劃分網格后就可以進行求解，此外該軟件具有強大的后處理功能，能夠對結果數據進行各種形式的處理并繪制圖像，便于研究人員對結果的分析[1]。

煤層氣俗稱瓦斯，是煤炭開采與含煤巖地區隧道工程的重大災害隱患，但又是一種潔凈的能源，我國煤層氣資源豐富,約，但勘探試驗工作起步較晚，煤層氣抽放過程中其運移過程與機理復雜，因此，深究煤層氣運移機理，對瓦斯抽采地面井的流場進行分析，合理開發煤層氣資源，不僅能為國民經濟的發展提供優質能源，有利于保護生態環境，而且有助于分析和了解瓦斯抽放的機理,尋求鉆孔或巷道抽放瓦斯的合理布置方式[2]，從根本上消除或減少瓦斯爆炸事故，這都具有十分重要的意義[3]。

在瓦斯抽放過程中，鉆孔周圍的瓦斯在多孔介質中的流動速度由慢速流變為快速流變，對于這種過渡流動的模擬，屬于需要自編程序計算的領域，因為轉換不同的流動定律就需要轉換不同的數學表達式。Darcy定律描述多孔介質中與管道有一定距離的慢速流動，Navier-Stokes方程應用于自由流動或明渠流動；在這兩者之間，存在多孔介質流，而其內剪切力不可忽略，就需要應用Brinkman或Forcheimer方程。本論文展示了如何應用COMSOL Multi-physic中的地球科學模塊內建的方程組來模擬這一流態轉換過程[4][5]。

1 模型的建立[6][7]

本算例模擬了煤層氣流向并流入一個鉆孔的流態轉換過程。首先，用耦合Darcy定律與Brinkman方程的方法，分析了多孔介質中的流動及其加速進入井筒的過程。然后，用耦合Darcy-Brinkman模型與Navier-Stokes方程的方法，模擬了進入井筒和其后的流體運動。與直覺相反，我們知道瞬態的Brinkman方程和Navier-Stokes方程相對來說比較容易求解。而本例采取的是穩態系統求解分析。

1.1 基本假設

圖1是Darcy-Brinkman模型的一個大致描述。煤層氣通過一個很薄的多孔介質層流，從一個橫向孔井流向井筒。在遠場(1m＜x＜4m)，流動遵循Darcy定律，而在井的側孔開口處（0.1m＜x＜1m），則遵循Brinkman方程。含氣層是0.875m厚，上下均為不可穿透的隔氣層，約束著有滲透性的含氣層。為簡化問題，假定含氣層具有各向同性的流體力學性質，并且流體的密度和粘性都是常量。進口處的流量和孔井處的壓強是已知的。流動為穩態的。

圖1 Darcy-Brinkman模型

1.2 流場控制方程

1.2.1 Darcy定律[8]

Darcy定律描述由壓力梯度和重力勢差引起的流動。Darcy定律的因變量是壓強,P。流動是緩慢的，因此動壓頭可以忽略。對于穩態情況，控制方程如下：

在這個方程中,K為滲透系數 （m2），η為動力學粘性系數kg/（m·s），ρf為流體密度（kg/m3），g為重力加速度（m/s2）。D表示垂直向上方向的坐標（m），Qs為相對于單位體積流體源儲量的體積流量（1/s）。在本問題中由于流場本身是一薄層，重力勢差可以忽略，因此D設為0。由于這一模型涉及多個流動定律，所以給壓力p加了下標dl，表示是Darcy定律的變量。

對于穩態分析來說，進入集水區的流量等于泵流量。Darcy定律的進口條件定義為：

式中：W為流向孔井方向的泵流量（m3/s）；為儲氣區半徑（m）；b為儲氣區厚度（m）。在Darcy流動和Brinkman流動之間的界面上，為了求得一個連續的解，Darcy定律和Brinkman方程的壓強和速度必須一致。因為在進口處已經用Neumann方式定義了通量，所以在Darcy-Brinkman界面上定義壓強：

其中下標“br”表示Brinkman方程。這個表達式給出了Dirichlet邊條件。

只有含氣層是可滲透到，上下均是不可滲透的隔氣層，Darcy定律的邊界條件定義為：

其中n是邊界上的法向單位矢量。

1.2.2 Brinkman方程[9]

Brinkman方程描述的是多孔介質中的高速流動，由于速度足夠高，因此內剪切力能夠引起動量傳播。Brinkman問題是x方向及y方向的動量平衡方程和連續性方程的聯立，其因變量有速度矢量分量u和v，以及壓強p。穩態的Brinkman方程為：

式中：ρ為密度（kg/m3）；η為動力學粘性系數（kg/（m·s））；u為速度矢量（m/s）；p為壓強（kg/（m·s））；ε為孔隙率,k（m2）為滲透系數。體積力項F（N/m3）可以包含微小的重力影響和壓縮效應，在本文中等于0。有文獻認為brinkman方程中的滲透系數k和Darcy定律中的滲透系數K是不同的,但是本文在兩個計算區域中應用同樣的滲透系數。

在Brinkman一側的Darcy-Brinkman界面上的約束條件是速度，因為在Darcy一側的約束條件是壓強。在界面上Brinkman一側的速度約束條件顯示了在Brinkman方程中速度是一個變量，這與Darcy定律不同。速度邊界條件寫為：

COMSOL Multi-physic可以自動地從Darcy方程計算出速度，將之應用于Brinkman方程的邊界條件。隔氣層和井壁對流體來說是不可滲透的。將其近似的定義為無滑移邊條件：

這一方程表示在邊界上一切速度分量都為0。

求此問題的特解要求在孔井處定義壓強，因為在其他邊界上，給出的都是通量邊條件。這個壓強的約束可以簡單地寫為：

對于Brinkman問題，邊界條件為

ubr=udlDarcy-Brinkman邊界條件

ubr=0 限制層邊界條件

ubr=0 井筒邊界條件

pbr=pwell穿孔邊界條件

在本文的分析中，兩個方程系統均是分析流體流動的壓強分布和速度，在根本上是可以相容的。但是，在Darcy定律中，因變量僅為壓強，而在Brinkman方程組中壓強和速度矢量均為因變量。因變量數量的不同造成了方程形式的輕微的不相容，可以用在方程系統中加入非理想弱項約束。弱項約束在方程系統中加入了新的積分方程，其拉格朗日乘子為因變量。弱項約束在Darcy定律內增加了一個拉格朗日乘子，在Brinkman方程組加入兩個，來彌補兩個方程系統在邊界上自由度數量的不同。增加拉格朗日乘子的步驟很簡單：打開應用模型屬性對話框，在“弱項約束”選單中選擇“開”。在“約束類型”中選擇“非理想”。

1.2.3 物性參數

表1是數值模擬物性參數取值情況。

表1 物性參數取值情況

2 建模及其求解

圖2為Darcy-Brinkman問題的模擬解，從遠場到孔井的緩慢流動過程中，應用的是Darcy定律；而孔井附近應用的則是Brinkman方程。兩種流動之間的轉換發生在x=1m處。流線表現了從右端的進口到左端的孔井的煤層氣流動。流線呈漏斗狀匯聚是因為氣體流向了井壁上的橫向孔洞。

圖3和圖4分別顯示了穿過Darcy-Brinkman界面，直到到2m距離處的壓強和速度的估算值。這兩者在Brinkman-Darcy界面上平滑過渡。隨著與井孔的距離增加，壓強升高并將流體推向井孔。而速度則隨著距離增加而減小，直至在Darcy流動區域內達到一個定值。由于從 Brinkman流動 (r＜1m)到Darcy(r＞1m)流動之間的壓強過渡平穩，因此可以求得一個穩態解。

圖2 數值模擬結果

圖3 壓力變化

圖4 速度變化

本論文展示了將兩種雖然有不同的因變量，但是可以相容的流動定律，Darcy定律和Brinkman流動方程聯立求解的簡單易懂的方法。這兩個流動定律描述的都是多孔介質流。因為Darcy定律沒有內剪切力引起的動量傳播項，在快速流動的區域內它就會把流量估計得過高。耦合Brinkman方程就可以描述這一額外的能量傳輸過程。

3 結 論

利用COMSOL Multi-physic軟件對煤層氣的流場進行數值模擬,模擬壓強和速度隨時間的變化結果如上圖3和圖4所示,從圖中可以看出,模擬數據與文獻[10]提供的實測數據具有較好的吻合特性，說明采用COMSOL Multi-physic軟件對所建立的煤層氣流場數學模型能夠很好的模擬。且與實際測試結果相符[10]；因此，通過壓力和速度曲線，可以分析得出煤層中壓力變化和速度變化的規律，從而可以為煤層甲烷的開采和預防瓦斯突出提供一定的參考。