基于COMSOL Multiphysics的十字狀裂隙中流體運動特征分析

樊嘉偉， 李紅星

(東華理工大學地球物理與測控技術學院，南昌 330000)

地震波在傳播過程中的耗散和衰減，一直以來都是很多工程領域所關心的問題，在以彈性波傳播為主要手段的勘探中，波的頻散與衰減十分普遍，其中頻散是指波在介質中的傳播速度會隨頻率發生改變，衰減則指的是波在介質中傳播時，能量被吸收導致的振幅減小。

彈性波勘探對裂隙性油氣儲層具有較強的敏感性，尤其波場的衰減及速度的頻散對于揭示地下儲層的性質具有良好的指示作用，因此彈性波作為油氣勘探的主要手段具有廣闊的發展前景[1-4]。地震波在地下介質中傳播時，眾多因素綜合造成了地震波的強衰減及高頻散現象，根據不同壓力梯度特征導致的流體流動形式可在空間尺度上將其分為3大類，分別為宏觀全局流動，又稱Biot全局流動，微觀噴射流動，又稱Squirt流動，和介觀尺度下流體的流動。地震波在含軟、硬孔隙斑塊飽和介質傳播過程中會誘發多個尺度孔隙流體流動而產生衰減和速度頻散,并且多個尺度間的流體流動相互影響[5-6]。很多學者已經提出大量的巖石物理理論，以此希望給出確切的理論依據來描述不同情況地震波的傳播機制。學界普遍認為造成微觀或介觀尺度下地震波能量衰減及速度頻散的主要原因是由于巖石骨架和孔隙流體的相對運動[1,7-12]。然而，目前還沒有能夠根據微觀或介觀模型的詳細幾何性質來解釋頻率相關衰減的實驗室測量結果的相關研究。

簡單噴射流的模型基于交互連接的水平和垂直扁平孔在地震波的激勵下，壓力梯度會引起垂直于波傳播方向的扁平裂隙中流體流動到臨近的平行于波傳播方向的裂隙中[13]。波動引起的流體在微觀尺度上的這類流動通常被稱為噴射流[13-18]。扁平毛孔在實際巖石中的表現形式通常為微裂紋和松散顆粒的互相接觸間隙。Gueguen等[19]根據微觀裂隙類型將其等效為3種不同的模型，即周期性和隨機分布的平面裂隙模型，以及硬幣形裂隙模型，并分析了不同模型地震波的振幅衰減及頻散情況。地震波激勵下介質中所發生的噴射流流動會隨頻率變化，在某一固定頻率時衰減值達到最大，即所謂的特征頻率，特征頻率受流體黏度和水平裂隙的幾何形狀控制。對于某些低滲性巖石，毛孔幾何形狀又與有效滲透率[20-21]密切相關。地震波激勵下孔隙中發生的噴射流作用對波的衰減具有不可忽略的影響，其衰減幅度可以在理想的幾何形狀孔隙作出很好的觀測[18]。Mavko等[14]提出單個十字狀裂隙的微觀噴射流機制，闡述了在地震波的激勵下，扁平裂隙中的流體會噴射到垂直裂隙中，他們的研究工作表明，當巖石中存在微觀裂隙時，即使含有極少量的流體，也會引起地震波強烈的衰減。聲波頻段，理想多孔介質中噴射流對波的衰減作用已經給出了良好的解釋[18]。Rubino等[21]研究表明，與噴射流模型中的扁平孔隙中的流動類似，在介觀尺度下，地震波大部分能量損失也發生在裂縫內，因為裂縫比巖石骨架背景相具有更高的滲透率，相關的衰減峰值與微觀裂隙的特征頻率相比更高[22-27]。隨后，Rubino等[22]又通過應用多孔彈性Biot方程對介觀尺度下噴射流機制進行數值模擬，該模型通過嵌入到低滲透性巖石骨架背景相的高滲裂隙來表示介觀尺度下實際巖石中噴射流發生的場所，他們認為在波的激勵下，巖石固相及液相的相對運動皆遵循Biot方程。很顯然，這種假設并不合理，因為Biot理論[28-30]的假設條件認為介質在宏觀上是統計各向同性的，介質中任意一點的變化都可用開放邊界的特征單元體來求取。而Rubino等的假設模型是固液相分離，且在不同代表性位置處固相滲透率并不一致，因此，在應用Biot方程對噴射流進行數值模擬時的實際意義并不嚴格。此外，由于缺乏切實有效的用于模擬巖石噴射流的數值模擬技術和微觀三維圖像中空間分辨率不足等問題，在實際巖石地震波引發裂隙結構中流體流動的研究工作中，沒有一種較優的方法來描述固液相的相互耦合機制。

基于上述缺陷，提出一種關于處理噴射流機制的新的數值方法，假設模型還是應用固液相分離的十字狀裂隙模型，噴射流發生場所在外部為極低滲透性固相，內部為流體飽和的十字狀交互裂隙內，并假設液相為不可壓縮流體。認為在地震波激勵下對于固相骨架部分的形變遵循應力應變平衡方程，流體在壓力作用下的流動遵循動量守恒的Navier-Stoke方程，固液兩相的相互耦合作用致使了地震波能量的衰減。并借助COMSOL Multiphysics軟件對地震波激勵下十字狀裂隙噴射流模型進行了數值模擬研究，通過對噴射流的發生過程及流體速度及壓力分布情況的分析，幫助學者們更好地觀察和理解地震波將能量傳遞給孔隙介質引發流體運動以及波能量衰減的動態過程。

1 理論依據

微觀噴射流的發生過程存在多個物理場之間的耦合，當地震波在巖石中傳播時，引起巖石結構變形，同時作為載荷將能量作用與孔隙內流體，從而引起流體的運動。巖石在微小外力作用下的應力變化遵循基本的平衡方程。

(1)

式(1)中：σ為應力張量；f為外加荷載，N。根據廣義胡克定律，對于各向同性材料，其應力張量分量滿足應力應變關系，即

σij=λεkkδij+2Gεij

(2)

對于彈性理論問題，以位移作為基本變量，在方程中消去應變和應力，可以得到用位移u(單位為m)表示的平衡方程，即Lame-Navier方程：

(λ+G)·u+G·u+f=0

(3)

孔隙流體的位移和速度分布遵循動量守恒的Navier-Stokes方程，質量守恒的連續性方程為

(4)

(5)

式中：p為流體應力張量；l為單位張量；μ為膨脹黏性系數，m2/s，一般情況下μ=0；s為變形速率張量，其在直角坐標中的分量為

(6)

2 數值模擬

利用COMSOL Multiphysics固體力學流固耦合模塊對地震波激勵下十字狀裂隙中流體流動形式進行數值模擬，COMSOL Multiphysics是瑞典COMSOL公司于1998年首次發布的一款基于有限元法的數值模擬軟件，前身是Matlab的一個工具箱，被稱為Toolbox 1.0。后來改名為Femlab 1.0(FEM為有限元，LAB是取用的Matlab)，經過多年的更新發展，COMSOL Multiphysics已成長為適用于仿真模擬工程、制造和科研等各個領域的研究型軟件。該軟件包含結構力學、流體流動、熱傳導、聲學、電磁學等模塊，具有強大的求解和后處理功能。COMSOL Multiphysics以定義和耦合任意數量偏微分方程的能力使之成為一個強大的分析工具，其靈活性和基于方程建模的形式可以幫助用戶進行更深入的研究。尤其是自定義偏微分方程(partial differential equation, PDE)應用模式，允許用戶自己建立模型，定義模型不同域所遵循的物理方程，通過求解瞬態或頻域下方程的解來達到數值分析的目的。在工程領域，COMSOL Multiphysics常被用于模擬地下巖層石油、天然氣資源、地下煤層和瓦斯資源、地下水資源開采以及污染物傳遞運移等。COMSOL Multiphysics軟件不僅可以實現多個物理場耦合問題的求解，同時還可以為用戶提供一個二次開發的優良環境。經過多年的開發，該軟件的模塊已經非常健全，對于大部分物理現象都有良好的仿真過程，同時，該軟件涵蓋了多種計算方法，對求解不同問題時，給用戶提供了更多的選擇，其計算精度高、速度快、分析功能強大、學科覆蓋范圍廣、操作簡單等特點贏得了廣大工程學者的青睞。對于筆者研究的問題，還可采用一些其他的正演模擬軟件進行分析，如美國GMS軟件中的MODFLOW模塊，德國FEFLOW軟件等，但目前國內還沒有較好的功能強大的數值模擬軟件可以進行詳細的分析。

在應用COMSOL Multiphysics進行模擬仿真時，主要應用其固體力學模塊中的流固耦合應用模式，該模塊經過多年發展和健全，已經擁有了完善的分析能力，針對筆者研究的問題，通過模型的建立和微分方程及參數設置，最后取得了良好的計算結果。數值模擬的主要步驟包括幾何模型的建立及材料選擇、網格剖分、激勵源的加載及邊界條件的設置、求解以及后處理等過程。以下將詳細介紹應用COMSOL Multiphysics對十字狀裂隙噴射流數值模擬工作流程，采用十字狀裂隙模型等效實際地層裂隙，如圖1所示。

圖1 實際地層等效模型Fig.1 Equivalent model of fractures

2.1 幾何模型的建立

由于COMSOL Multiphysics有很多外部接口，對于大部分格式都有良好的兼容性，因此在建立模型時有很多方法，可以直接應用軟件中固有的幾何開發工具欄來創立模型，設置好模型參數，也可以通過CAD、MAPGIS等其他繪圖軟件提前構建好模型，再導入到COMSOL Multiphysics中，對用戶提供了多種選擇，因此應用起來十分方便。

由于所采用的十字狀裂隙模型相對簡單、規則，因此應用軟件自帶幾何工具欄進行模型的建立。操作步驟如下：

(1)打開COMSOL Multiphysics軟件，點擊二維，進入物理場選擇界面，然后單擊流體流動欄，選定流固耦合(fsi)模塊，研究模式設置為瞬態，然后點擊完成。

(2)在全局菜單下定義、選擇、輸入或導入模型參數以及理論解析公式參數，并給出相應單位。

(3)在幾何菜單下開始建立十字狀裂隙模型，如圖2所示。

(4)建立好幾何模型后，開始對不同的域定義材料。添加材料的過程有兩種，可以直接在COMSOL Multiphysics內置材料庫中直接選擇相應的材料，也可以通過添加空材料，根據所求的物理場，直接設置需要的材料屬性參數，直接從內置材料庫中直接選取材料，其固體域和液體域分別為花崗巖和水。材料屬性如圖3、圖4所示。

圖2 十字狀裂隙模型Fig.2 Model of interconnectedfractures

圖3 固體域材料屬性Fig.3 Material properties of the solid domain

圖4 流體域材料屬性Fig.4 Material properties of the liquid domain

2.2 網格剖分

對模型的剖分如圖5所示，應該注意的是，所添加的激勵源是一個高頻短時脈沖。為了更好地刻畫該脈沖信號在介質中的傳播，保證其有足夠多的網格節點，因此在模型網格化時，應根據脈沖信號頻率和材料屬性計算出最小網格單元尺寸。

圖5 十字狀裂隙模型的網格剖分Fig.5 Meshing of the interconnected fractures model

同時，用戶自定義網格時還應該注意對于二維模型通常采用三角形或四邊形網格剖分，三維模型則采用四面體或六面體網格剖分。除此之外，應該注意模型中不同位置處網格的疏密程度，如在流固耦合邊界處應適當加密網格，有利于提高計算精度，使計算結果更符合實際情況，在其他位置則可以適當稀疏網格，不僅可以使得計算收斂性好，而且能夠提升計算速度。

2.3 激勵源的加載及邊界條件的設置

在應用COMSOL Mutltiphysics進行模擬仿真時，通過添加一個短時脈沖的載荷作為激勵源來模擬地震波在介質中的傳播過程。由于模型相對較小，可以認為地震波是以平面波的形式入射并通過該模型的。因此，激勵源的加載是在其上界線位置處添加一個邊界載荷來實現。在彈性波激勵的外力作用下，為了更符合實際情況中地震波在地下巖石中的傳播過程，期望模型不會在垂直方向上發生位移，因此在下界面位置處，應該添加一個固定約束作為條件。兩側邊界則可認為是自由邊界。流體和固體界面無需設置，系統內部會自動默認為流固耦合界面。圖6所示為本次實驗所采用的激勵源，其表達式為r(t)=A0[1-2(ωπt)2]e-(ωπt)2的地震子波，ω為主頻，本次實驗主頻取30 kHz，A0為子波的振幅，t為時間。

圖6 主頻為30 kHz的激勵源Fig.6 Excitation source with a frequency` of 30 kHz

2.4 求解及后處理

2.4.1 求解

在求解過程中，由于地震波的傳播途徑是貫穿整個模型的，因此應將最小研究時間設置為波從模型上界面傳播到下界面的走時，時間步長可以根據實際需要定義。選取的研究時長和步長如圖7所示。

圖7 研究時間和步長的設定Fig.7 Settings of study duration and stepping

2.4.2 后處理

最后可以將數值模擬所得到的計算結果根據研究需要，分別繪制成相應的一維、二維圖件。COMSOL Multiphysics軟件自身具有非常全面和強大的繪圖功能，繪圖之前，需要在數據集中設置好要觀測的對象，例如繪制某一時刻某一點處的任意物理場結果，任意點處任意時刻的物理場結果等，根據需要，繪制不同時刻在水平和垂直裂隙中液體的流動速度場和壓力場。

3 結果分析

圖8所示為主頻為30 kHz時地震波在含流體孔隙介質中傳播時不同時刻的波場快照，從圖8中可以看出，地震波是以平面波的形式入射的，在介質中傳播時受到裂隙壁面影響發生了一定程度的反射。

圖8 地震波在不同時刻介質中的波場快照Fig.8 Wave field snapshot of seismic waves at different times

圖9所示為地震波作用下，不同時刻十字狀裂隙內流體的壓力及速度分布情況，x和y代表的坐標軸皆表示預設模型水平及垂直刻度。圖9中可以直觀看到顏色越深的部分場值越大。

圖9 不同時刻流體的速度及壓力場Fig.9 Velocity and pressure field of fluid at different times

圖10 水平方向選取的二維截線Fig.10 Cut line in horizontal fracture

圖11 垂直方向選取的二維截線Fig.11 Cut line in vertical fracture

為了觀測水平裂隙及垂直裂隙中流體的速度及壓力分布情況，在水平裂隙和垂直裂隙的軸向上分別設置了一條截線，如圖10、圖11所示。得到不同情況下的流體速度及壓力分布情況如圖12～圖15所示。

圖12 水平截線上不同時刻流體速度Fig.12 Fluid velocity at different times on the horizontal line

圖13 水平截線上不同時刻流體壓力Fig.13 Fluid pressure at different times on the horizontal line

圖14 垂直截線上不同時刻流體速度Fig.14 Fluid velocity at different times on the vertical line

圖15 垂直截線上不同時刻流體壓力Fig.15 Fluid pressure at different times on the vertical line

由圖12、圖13可以看到，隨著時間變化，在水平裂隙中，不同時刻，越靠近中心位置，流體的流速越大，壓力越小，極值都在裂隙的中心位置處；由圖14可以看到，不同時刻在垂直裂隙中越靠近中心位置，流速也是逐漸增大，在中心位置處流體速度達最大值；在圖15中可以看到，不同時刻，壓力極值也是出現在中心位置處，但由于壓縮波的傳播方向與垂直裂隙的軸向一致，在波傳播方向上使流體強制運動，另一方面，水平裂隙中流體在壁面壓力下流入到垂直裂隙中，造成了流體壓力不規則變化的現象。

4 結論

應用COMSOL Multiphysics仿真軟件對十字狀裂隙噴射流機制進行了數值模擬，得到了在地震波激勵下，不同時刻水平及垂直裂隙中流體的壓力及速度分布情況。在水平裂隙中，越靠近中心位置處，流體流速越大，壓力越小，在中心處存在極值；在垂直裂隙中，靠近中心位置，流體流速也是逐漸增大，在中心位置處存在極值，但由于波傳播方向和水平裂隙流體擠壓的影響，導致垂直裂隙中流體壓力隨時間出現不規則的變化。通過分析和對比可知，所得結果符合噴射流機制中流體的運動規律，與前人相關研究比較，雖然存在差異，但大致情況基本吻合。

此外，本次實驗對該物理問題進行仿真時，應用COMSOL Multiphysics仿真軟件取得了較好的結果，證明COMSOL Multiphysics對該問題研究的可靠性。地震波引發介質中流體流動導致的耗散問題一直以來都是工程界所關心的問題，前人對此做了大量的研究工作，也提出了很多適用于不同情況下的理論，但直至目前為止，還沒有一個能夠適用于所用情況的統一的理論。在自然條件下，地下介質分布情況非常復雜，不能將地震波的頻散或衰減簡單地認為是由某一個或一些因素決定的，因此很難通過所觀測的數據針對某一特定問題進行單獨研究。然而借助COMSOL Multiphysics軟件，通過設置不同的模型，可以對所關心的內容進行單獨分析，其強大的功能能夠直觀便捷地反映所研究的問題，為科研人員提供了一個良好的工作平臺。