翟清偉
(兗州煤業(yè)股份有限公司鮑店煤礦,山東 濟寧 273500)
COMSOL Multi-physic是一種多重物理量耦合軟件,該軟件是由MATLAB軟件工具箱發(fā)展而來的,以有限元方法進行分析求解,其優(yōu)點在于高度的靈活性,強大的求解器和較高的計算精度,進行求解時只需要將所建立的數(shù)學模型輸入軟件的PDE模塊中,設定求解域,指定邊界條件并劃分網(wǎng)格后就可以進行求解,此外該軟件具有強大的后處理功能,能夠?qū)Y(jié)果數(shù)據(jù)進行各種形式的處理并繪制圖像,便于研究人員對結(jié)果的分析[1]。
煤層氣俗稱瓦斯,是煤炭開采與含煤巖地區(qū)隧道工程的重大災害隱患,但又是一種潔凈的能源,我國煤層氣資源豐富,約,但勘探試驗工作起步較晚,煤層氣抽放過程中其運移過程與機理復雜,因此,深究煤層氣運移機理,對瓦斯抽采地面井的流場進行分析,合理開發(fā)煤層氣資源,不僅能為國民經(jīng)濟的發(fā)展提供優(yōu)質(zhì)能源,有利于保護生態(tài)環(huán)境,而且有助于分析和了解瓦斯抽放的機理,尋求鉆孔或巷道抽放瓦斯的合理布置方式[2],從根本上消除或減少瓦斯爆炸事故,這都具有十分重要的意義[3]。
在瓦斯抽放過程中,鉆孔周圍的瓦斯在多孔介質(zhì)中的流動速度由慢速流變?yōu)榭焖倭髯?,對于這種過渡流動的模擬,屬于需要自編程序計算的領域,因為轉(zhuǎn)換不同的流動定律就需要轉(zhuǎn)換不同的數(shù)學表達式。Darcy定律描述多孔介質(zhì)中與管道有一定距離的慢速流動,Navier-Stokes方程應用于自由流動或明渠流動;在這兩者之間,存在多孔介質(zhì)流,而其內(nèi)剪切力不可忽略,就需要應用Brinkman或Forcheimer方程。本論文展示了如何應用COMSOL Multi-physic中的地球科學模塊內(nèi)建的方程組來模擬這一流態(tài)轉(zhuǎn)換過程[4][5]。
本算例模擬了煤層氣流向并流入一個鉆孔的流態(tài)轉(zhuǎn)換過程。首先,用耦合Darcy定律與Brinkman方程的方法,分析了多孔介質(zhì)中的流動及其加速進入井筒的過程。然后,用耦合Darcy-Brinkman模型與Navier-Stokes方程的方法,模擬了進入井筒和其后的流體運動。與直覺相反,我們知道瞬態(tài)的Brinkman方程和Navier-Stokes方程相對來說比較容易求解。而本例采取的是穩(wěn)態(tài)系統(tǒng)求解分析。
圖1是Darcy-Brinkman模型的一個大致描述。煤層氣通過一個很薄的多孔介質(zhì)層流,從一個橫向孔井流向井筒。在遠場(1m<x<4m),流動遵循Darcy定律,而在井的側(cè)孔開口處(0.1m<x<1m),則遵循Brinkman方程。含氣層是0.875m厚,上下均為不可穿透的隔氣層,約束著有滲透性的含氣層。為簡化問題,假定含氣層具有各向同性的流體力學性質(zhì),并且流體的密度和粘性都是常量。進口處的流量和孔井處的壓強是已知的。流動為穩(wěn)態(tài)的。
圖1 Darcy-Brinkman模型
1.2.1 Darcy定律[8]
Darcy定律描述由壓力梯度和重力勢差引起的流動。Darcy定律的因變量是壓強,P。流動是緩慢的,因此動壓頭可以忽略。對于穩(wěn)態(tài)情況,控制方程如下:
在這個方程中,K為滲透系數(shù) (m2),η為動力學粘性系數(shù)kg/(m·s),ρf為流體密度(kg/m3),g為重力加速度(m/s2)。D表示垂直向上方向的坐標(m),Qs為相對于單位體積流體源儲量的體積流量(1/s)。在本問題中由于流場本身是一薄層,重力勢差可以忽略,因此D設為0。由于這一模型涉及多個流動定律,所以給壓力p加了下標dl,表示是Darcy定律的變量。
對于穩(wěn)態(tài)分析來說,進入集水區(qū)的流量等于泵流量。Darcy定律的進口條件定義為:
式中:W為流向孔井方向的泵流量(m3/s);為儲氣區(qū)半徑(m);b為儲氣區(qū)厚度(m)。在Darcy流動和Brinkman流動之間的界面上,為了求得一個連續(xù)的解,Darcy定律和Brinkman方程的壓強和速度必須一致。因為在進口處已經(jīng)用Neumann方式定義了通量,所以在Darcy-Brinkman界面上定義壓強:
其中下標“br”表示Brinkman方程。這個表達式給出了Dirichlet邊條件。
只有含氣層是可滲透到,上下均是不可滲透的隔氣層,Darcy定律的邊界條件定義為:
其中n是邊界上的法向單位矢量。
1.2.2 Brinkman方程[9]
Brinkman方程描述的是多孔介質(zhì)中的高速流動,由于速度足夠高,因此內(nèi)剪切力能夠引起動量傳播。Brinkman問題是x方向及y方向的動量平衡方程和連續(xù)性方程的聯(lián)立,其因變量有速度矢量分量u和v,以及壓強p。穩(wěn)態(tài)的Brinkman方程為:
式中:ρ為密度(kg/m3);η為動力學粘性系數(shù)(kg/(m·s));u為速度矢量(m/s);p為壓強(kg/(m·s));ε為孔隙率,k(m2)為滲透系數(shù)。體積力項F(N/m3)可以包含微小的重力影響和壓縮效應,在本文中等于0。有文獻認為brinkman方程中的滲透系數(shù)k和Darcy定律中的滲透系數(shù)K是不同的,但是本文在兩個計算區(qū)域中應用同樣的滲透系數(shù)。
在Brinkman一側(cè)的Darcy-Brinkman界面上的約束條件是速度,因為在Darcy一側(cè)的約束條件是壓強。在界面上Brinkman一側(cè)的速度約束條件顯示了在Brinkman方程中速度是一個變量,這與Darcy定律不同。速度邊界條件寫為:
COMSOL Multi-physic可以自動地從Darcy方程計算出速度,將之應用于Brinkman方程的邊界條件。隔氣層和井壁對流體來說是不可滲透的。將其近似的定義為無滑移邊條件:
這一方程表示在邊界上一切速度分量都為0。
求此問題的特解要求在孔井處定義壓強,因為在其他邊界上,給出的都是通量邊條件。這個壓強的約束可以簡單地寫為:
對于Brinkman問題,邊界條件為
ubr=udlDarcy-Brinkman邊界條件
ubr=0 限制層邊界條件
ubr=0 井筒邊界條件
pbr=pwell穿孔邊界條件
在本文的分析中,兩個方程系統(tǒng)均是分析流體流動的壓強分布和速度,在根本上是可以相容的。但是,在Darcy定律中,因變量僅為壓強,而在Brinkman方程組中壓強和速度矢量均為因變量。因變量數(shù)量的不同造成了方程形式的輕微的不相容,可以用在方程系統(tǒng)中加入非理想弱項約束。弱項約束在方程系統(tǒng)中加入了新的積分方程,其拉格朗日乘子為因變量。弱項約束在Darcy定律內(nèi)增加了一個拉格朗日乘子,在Brinkman方程組加入兩個,來彌補兩個方程系統(tǒng)在邊界上自由度數(shù)量的不同。增加拉格朗日乘子的步驟很簡單:打開應用模型屬性對話框,在“弱項約束”選單中選擇“開”。在“約束類型”中選擇“非理想”。
1.2.3 物性參數(shù)
表1是數(shù)值模擬物性參數(shù)取值情況。
表1 物性參數(shù)取值情況
圖2為Darcy-Brinkman問題的模擬解,從遠場到孔井的緩慢流動過程中,應用的是Darcy定律;而孔井附近應用的則是Brinkman方程。兩種流動之間的轉(zhuǎn)換發(fā)生在x=1m處。流線表現(xiàn)了從右端的進口到左端的孔井的煤層氣流動。流線呈漏斗狀匯聚是因為氣體流向了井壁上的橫向孔洞。
圖3和圖4分別顯示了穿過Darcy-Brinkman界面,直到到2m距離處的壓強和速度的估算值。這兩者在Brinkman-Darcy界面上平滑過渡。隨著與井孔的距離增加,壓強升高并將流體推向井孔。而速度則隨著距離增加而減小,直至在Darcy流動區(qū)域內(nèi)達到一個定值。由于從 Brinkman流動 (r<1m)到Darcy(r>1m)流動之間的壓強過渡平穩(wěn),因此可以求得一個穩(wěn)態(tài)解。
圖2 數(shù)值模擬結(jié)果
圖3 壓力變化
圖4 速度變化
本論文展示了將兩種雖然有不同的因變量,但是可以相容的流動定律,Darcy定律和Brinkman流動方程聯(lián)立求解的簡單易懂的方法。這兩個流動定律描述的都是多孔介質(zhì)流。因為Darcy定律沒有內(nèi)剪切力引起的動量傳播項,在快速流動的區(qū)域內(nèi)它就會把流量估計得過高。耦合Brinkman方程就可以描述這一額外的能量傳輸過程。
利用COMSOL Multi-physic軟件對煤層氣的流場進行數(shù)值模擬,模擬壓強和速度隨時間的變化結(jié)果如上圖3和圖4所示,從圖中可以看出,模擬數(shù)據(jù)與文獻[10]提供的實測數(shù)據(jù)具有較好的吻合特性,說明采用COMSOL Multi-physic軟件對所建立的煤層氣流場數(shù)學模型能夠很好的模擬。且與實際測試結(jié)果相符[10];因此,通過壓力和速度曲線,可以分析得出煤層中壓力變化和速度變化的規(guī)律,從而可以為煤層甲烷的開采和預防瓦斯突出提供一定的參考。