求解三維裂縫型多孔介質流體的場分裂預條件子研究

楊念 楊海建,* 邵柏強

(1.湖南大學重慶研究院,重慶,401120?2.湖南大學數學學院,長沙,410082)

1 引言

由于采用高分辨率網格進行流動模擬可以捕捉精細尺度現象,因此它常用于對復雜地質流體的物理現象進行精確描述.例如裂縫型油藏問題,在模擬過程中連通的裂縫只占據很小的計算區域,但可以完全控制從毫米到數百公里的流動過程[1].為了模擬長時間的移動過程和多尺度的復雜地形,往往需要利用非常多的網格塊進行全尺度計算[2,3],這使得完成預測可能需要數周甚至數月的時間[4,5,6].對于如此復雜、非線性的流動問題進行大規模模擬,需要的計算資源非常大.因此,為了提高計算效率必須建立能夠充分利用大規模并行計算機優勢的并行流場模擬器或求解器.

對于大型非線性方程組,Newton-Krylov 算法[7,8,9]是一種具有快速局部收斂特性的通用迭代求解方法.這種迭代方法可以在一些流行的開源庫中找到,如PETSc[10],NITSOL[11].由于Jacobi矩陣的高度非對稱和病態,在求解非線性系統時的一個主要瓶頸是如何在每次Newton 迭代時保證對線性系統預處理的魯棒性和高效性.如果沒有適當的預條件子,Newton-Krylov 算法可能不會有良好的收斂速度.因此,即使有超級計算機系統和大量的計算核心可用,也不能減少計算時間.一種重要的Newton-Krylov 預處理器被稱為場分裂(FS)預條件子,它將原始系統劃分為一系列簡化的隱式系統.FS 預條件子具有塊結構,使用物理變量作為拆分策略,然后應用不同的預處理器來處理具有不同屬性的變量.因此,用不同方法從原始系統矩陣中提取塊,可以得到不同的預條件子,如加性FS 預條件子[12,13]、乘性FS 預條件子[14]、舒爾補FS 預條件子[13,15]和約束壓力殘差(CPR)預條件子[16,17,18].加性或乘性FS 方法是基于Jacobi 或Gauss-Seidel 方法,分別將耦合問題分解為一系列更簡單的問題.舒爾補FS 方法也是為利用不同耦合過程的特性設計的.約束壓力殘差預條件子是Wallis 在1980 年代初提出的[17],是一個行業標準的求解器,已成功應用于各種流動問題[12,13,14,18,19],包括油藏模擬和計算流體動力學.

本文研究基于物理和區域分解方法的不同組合的FS 預條件子,應用于裂縫型多孔介質中的非定常流體問題[1,4].在該算法中,FS 預條件子中的所有子系統都采用基于RAS 的具有多層重疊的域分解技術進行求解,構造包括稀疏LU 或ILU 分解在內的不同子域求解器,并對其進行詳細比較.為了進一步提高算法的魯棒性和可擴展性,并節省計算時間,本文利用幾何多重網格框架構建一種兩水平FS 預條件子,通過粗網格校正,獲得粗網格和細網格之間附加的子域間信息,以解析誤差的低頻分量.數值實驗表明,所提出的預處理方法在超級計算機上對某些地下流動問題具有高度的可擴展性.

本文安排如下:在第2 節中介紹裂縫型多孔介質的雙孔雙滲模型?在第3 節中,給出基于物理和區域分解方法的不同組合的FS 預條件子,并將對應的算法應用于求解所得到的非線性系統?在第4 節中進行數值實驗來檢驗所建立算法的效率和有效性?最后在第5 節對本文進行歸納總結.

2 裂縫型多孔介質雙孔雙滲模型

本節介紹描述裂縫型多孔介質中的油藏流動的雙孔雙滲(DPDP)模型[1,4].在該模型中,基質和裂縫被視為兩個獨立的連續體,通過界面相互連接并允許它們轉移.DPDP 模型的質量平衡方程為

其中,?α為基質和裂縫的孔隙率,ρα為流體密度,uα為達西速度,Kα為滲透率張量,pα為多孔介質的壓力,μ為流體的粘度,g為重力矢量,計算域為??Rd,最終模擬時間為Te.此外,τ和q分別為轉移項和生產項,計算式如下:

其中,Lx,Ly和Lz分別表示x,y,z方向上的裂縫間距,θ是一個與生產井位置有關的常數,re是排水半徑,rw是井半徑,pb是井底壓力,pv是井周圍的平均裂縫壓力.

利用Peng-Robinson 狀態方程(PR-EoS,[5,20])將多孔介質中可壓縮流體ρα的密度描述為溫度和壓力的函數:

其中,W,R和T分別為分子量、通用氣體常數和溫度.氣體壓縮系數Zα由下式計算:

這里,PR-EoS 參數Aα和Bα定義如下:

其中,Tc和pc分別為氣體臨界狀態下的溫度和壓力,b是如下的分段函數:

其中,Tb是沸點溫度,patm是大氣壓.

針對上述問題,本文使用混合有限元(Mixed Finite Element,MFE)方法進行空間離散,用顯式第一步、單對角系數、對角隱式Runge-Kutta 格式進行時間離散(詳見參考文獻[21]).式(2.1)的完全離散形式是一個具有未知數(p1,p2)的非線性代數方程組:

3 Field-Split 預條件子

本節介紹求解非線性代數系統的Newton-Krylov 算法,并在Newton-Krylov 算法中,構建幾種基于物理和區域分解方法的的場分裂(FS)預條件子.

3.1 Newton-Krylov 算法

首先介紹求解非線性代數方程組(2.2)的Newton-Krylov 算法[7,8,9],記未知量組成的向量為X=(x1,x2,···,xl)∈Rn,殘差向量為F(X)=(F1(X),F2(X),···,Fi(X))T∈Rn.設(2.2)對應的Jacobi 矩陣為:

假設初始向量為X0∈Rn,第k次Newton 迭代的近似解為Xk,則第k+1 次Newton 迭代的近似解為:

其中,dk是用Krylov 迭代法[22]求得的搜索方向,λk∈(0,1]是由三次線搜索確定的步長.由此可得到求解(2.2)的迭代形式:

其中,?F（Xk）為Jacobi 矩陣(3.1),滿足以下條件:

其中,ηr∈[0,1)是線性迭代的相對誤差,ηa∈[0,1)是線性迭代的絕對誤差.迭代過程的終止條件為:

其中,εr和εa分別是Newton 迭代的相對誤差和絕對誤差.

在Newton-Krylov 算法中,求解器最重要的部分是選擇合適的預條件形式:

3.2 場分裂預條件子

場分裂預條件子的方法是使用域松弛或域分解來求解l×l線性塊系統.為簡單起見,將線性塊系統限制為2×2 塊,即(3.1)中的Jacobi 矩陣為:

本文給出四種場分裂預條件子,即加性FS 預條件子、乘性FS 預條件子、舒爾補FS 預條件子和CPR 預條件子.

加性FS 預條件子和乘性FS 預條件子分別為:

舒爾補FS 預條件子為:

其中,C=J22-.舒爾補FS 預條件子的逆可以寫成以下形式:

并且舒爾補FS 預條件子可以全分解為LDU形式,通過去掉某些項得到預條件子的幾種變體,如下三角(LD)-1,上三角(DU)-1和對角線等類型.注意,對角線“diag”類型具有以下形式:

在實際中,C的近似值有以下幾種形式:“a11”類型C-1=?“selfp” 類型C-1=(J22-J21(diag(J11))-1J12)-1?“full”類型C-1=.

第四種FS 預條件子是約束壓力殘差預條件子(CPR 預條件子).利用CPR 預條件子求解方程是一個兩階段過程:首先通過消除某些變量求解一個較小的簡化系統,然后近似求解原始的完整系統.第一階段預條件子為:

在兩水平FS 預條件子中,粗網格主要用于信息交換.在定義的粗網格上需要求解以下線性問題:

由于一些直接求解方法在大規模仿真應用中計算成本較大,所以采用一水平FS 預條件子的GMRES 方法.當用迭代法求解粗線性系統時,整體的預條件子是一個迭代過程,且預條件子在每次迭代中都有變化.因此,使用兩水平FS 預條件子求解線性系統時,采用柔性GMRES(fGMRES)方法.

3.3 限制加性的Schwarz 預條件子

在FS 預條件子中,用限制加性Schwarz(RAS)算法來近似對應變量Jacobi 矩陣的逆.假設待求解的原始系統AX=b有N個未知數和N個線性方程.在定義RAS 預條件子之前,先介紹區域? 的一種分解方式.首先將區域? 分解為Np個互不重疊的子區域?i,i=1,2,3,···,Np(Np為處理器的個數),然后將每個子區域?i向外擴張δ層網格,得到重疊子區域?i,δ,其中δ為子區域的重疊參數.圖1是區域? 的某一分解方式,其中灰色實線表示網格線,黑色實線表示區域? 的非重疊劃分方式,深灰色區域為某一非重疊子區域?i,0,淺灰色區域為該非重疊子區域?i,0向外擴張的δ=2 層網格.設Ni和N分別是區域和? 內的自由度.定義限制性算子Ri,δ∈,該矩陣的元素定義為:若1≤l1≤Ni和1≤l2≤N,對應于子區域?i,δ內的同一個網格點,則該元素的值取為1,否則為0.

限制加性Schwarz 預條件子定義為:

其中,Ri,0的定義與Ri,δ類似,當且僅當其對應的網格點在區域?i內時,矩陣元素的值為1.(Ai,δ)-1是子域Jacobi 矩陣Ai,δ的近似逆,它是由稀疏LU 或ILU 因子分解[22]計算出來的.

4 數值實驗

本節將進行一系列數值實驗來評估所提出的場分裂預條件子的性能.本文研究的算法使用開源的可移植、可擴展的科學計算工具包(PETSc)[10]來實現,數值測試在天河2 號超級計算機上進行.在以下測試的表格中,“N.It”表示每個時間步的Newton 迭代次數的平均值,“L.It”表示每次Newton 迭代一水平或兩水平FS 預處理GMRES 迭代的平均次數,“Time”是總計算時間(以秒為單位),“Np”表示處理器核心的數量.

4.1 算例一:拋物型滲透率

在算例一中,油藏的計算域是?=(0 m,50 m)3.如圖2(a)所示,在基體中,高滲透帶位于(y≥0.16x2-7.78x+112.22)[23],滲透率為100 md(圖中紅色部分),其他為低滲透帶,滲透率為1 md?裂縫中的滲透率值比基體中的滲透率值大100 倍,裂縫孔隙率為0.001,基體孔隙率為0.05?初始壓力隨深度的變化而變化,設定在10～11Mpa 之間?注水井和生產井的井底壓力分別為12Mpa和8Mpa.圖(b)–(d)為不同時刻裂縫壓力分布圖.

圖2 算例一中滲透率分布圖和不同時刻裂縫壓力分布圖

表1顯示了加性FS 預條件子和乘性FS 預條件子在不同子域求解器和重疊因子下的測試結果能.在測試中,網格大小固定為64×64×64,處理器的數量為Np=8,時間步長?t=0.1 天,計算5 個時間步.測試結果表明:(a)LU 子求解器在線性迭代方面優于ILU 子求解器,而ILU 子求解器在計算時間方面有優勢?(b)重疊因子越大,線性迭代次數越少,但由于通信成本增加,計算時間不一定減少.為線性迭代和計算時間之間的平衡,在接下來的測試中,我們固定一個最優參數組合為重疊因子δ=1 和求解器ILU(1).

表1 算例一中不同子域求解器和重疊因子對加性和乘性FS 預條件子的影響

表2顯示了在“full”類型的近似塊分解下,舒爾補FS 預條件子在不同預處理類型下的性能.而表3顯示了在舒爾補FS 的預處理類型固定為“a11”時,使用不同的近似塊分解時預條件子的性能的測試結果.測試結果表明:(a)線性迭代次數對網格尺寸不敏感?(b)預處理類型和近似塊分解的選擇對舒爾補FS 方法的魯棒性和效率有很大影響?(c)在舒爾補FS 預條件子中,“a11”預處理類型和“full”類型近似塊分解的組合是一個用于接下來的數值試驗的最優選擇.舒爾補FS 預條件子的線性迭代次數非常少,主要原因是“lower”類型或“upper”類型的舒爾補FS 預條件子具有1 或2次最小多項式,一般在2 次線性迭代后收斂.從圖4從可以看出,舒爾補FS 預條件子可以將Jacobi矩陣的特征值很好地限制在1 附近,因而線性迭代次數少.

表2 算例一中不同預處理類型對舒爾補FS 預條件子的影響(“?”表示計算兩小時仍無結果)

表3 算例一中不同近似塊分解對舒爾補FS 條件子的影響(“?”表示計算兩小時仍無結果)

4.2 算例二:階梯型滲透率

在算例二中,油藏的計算域是?=(0 m,800m)×(0 m,800m)×(0 m,400m).如圖3所示,基質中的紅色區域為高滲透區,滲透率為100md,藍色區域為低滲透區,滲透率為1md?裂縫中的滲透率值比基質中的滲透率值大100 倍?基質中的孔隙率值為0.2,比裂縫中的孔隙度大10 倍.注水井和生產井的井底壓力分別為12Mpa 和8Mpa.

圖3 算例二中滲透率分布圖和壓力分布圖

采用一水平的FS 預條件子求解上述DPDP 模型,表4顯示了幾種不同FS 預條件子在不同網格尺寸下的線性迭代和計算時間方面的性能.表中“CPR-p1”表示第一級預處理條件子(3.8)中的構造基于基質壓力p1,依此類推.從表4可以看出:(a)舒爾補FS 預條件子需要的GMRES 迭代次數最少?(b)加性FS、乘性FS 和CPR-p2預條件子的性能在線性迭代方面具有相似的數值結果?(c)平衡線性迭代和計算時間最優預條件子為CPR-p1.

表4 算例二中不同FS 預條件子性能的比較

圖4使用16×16×8 的網格繪制了該問題的原始Jacobi 矩陣和具有不同FS 預條件子的預處理矩陣的特征值.從圖中可以看出原始Jacobi 矩陣的特征值是非常分散的,在使用場分裂預條件子處理后特征值的分布聚集度高,并且CPR-p2預條件子的效果最好.

圖4 算例二中Jacobi 矩陣和預處理矩陣的特征值分布圖

4.3 算例三:SPE10 基準算例

算例三采用SPE10 基準算例(SPE10)[24]來測試三維裂縫型油藏模擬問題.由于滲透率和孔隙度的高度非均質性,該測試算例是油藏模擬中一個具有挑戰性的基準問題.如圖5所示,基質的滲透率在垂直方向的變化從6.65×10-4md 到2×104md 不等,其變化速度快于水平方向,基質孔隙度不再是一個常數,其范圍為0 到0.5.實驗中裂縫滲透率值比基質處大400 倍,裂縫孔隙度比基質小30 倍,基質和裂縫中的初始壓力隨深度變化,設定為10.0Mpa 到10.4Mpa.井底壓力為p=3.45Mpa.在數值試驗中采用60×220×85 網格對油藏區域1200ft×2200ft×170ft 進行模擬,從區域左側注入流量,從右側產出.

圖5 算例三中SPE10 基質滲透率對數分布圖和孔隙度分布圖

首先測試當時間步長?t改變時,所提出的不同FS 預條件子的性能,測試結果見表5.在測試中,模擬時間為T=1 天,式(3.3)中的線性停止條件為ηr=10-8和ηa=10-5,式(3.4)中的非線性停止條件為εr=10-10和εa=10-6.表5中的結果表明,采用不同FS 預條件子在幾種不同時間步長的迭代過程都是收斂的,且是無條件穩定的?在計算時間方面,CPR 預條件子的性能優于其他FS預條件子,并且隨著時間步長?t的減小,非線性和線性迭代的平均次數變小,而總計算時間增加,這與全隱式方法在各種應用中的結果一致.

表5 不同時間步長對FS 預條件子性能的影響

4.4 算例四:隨機滲透率

在算例四中,使用開源工具箱MRST 生成了網格大小為128×128×128 的隨機非均勻滲透率.如圖6所示,在基質中,隨機滲透率取對數后的分布范圍為[0.59,4.1],裂縫滲透率值比基質處大400 倍,基質中的孔隙度為0.05,裂縫中的孔隙度為0.001.在數值試驗中對油藏區域[0,100]3進行模擬,從區域左側注入流量,從右側產出.試驗固定模擬時間為T=1 天,數值結果見表5.從表中可以看出不同FS 預條件子在幾種不同時間步長的迭代過程都是收斂的.

圖6 算例四中隨機滲透率對數分布圖

4.5 并行測試

本小節測試所提出的并行求解器在將Newton-Krylov 算法與一水平或兩水平FS 預條件子相結合時的強可擴展性和弱可擴展性.在兩水平FS 預條件子中,粗細網格比固定為4,粗網格上的線性系統采用一水平RAS 預條件子的GMRES 方法求解,求解器絕對誤差為0.2.可擴展性的衡量指標是Speedup=T1/T2,其中T1和T2分別是使用Np,1和Np,2個處理器(Np,1≤Np,2)時,算法的計算時間.在進行強可擴展性的測試時,將時間步長固定為?t=0.1 天,網格設置為256×256×256,處理器的數量從128 個依次增加到2048 個,對算例一的前五個時間步長的情形進行模擬.表6為強可擴展性測試的非線性和線性迭代的次數以及在處理器數量下的計算時間,結果表明算法的總計算時間隨著處理器數量的增加而減少.在弱可擴展測試中,網格規模從128×128×128(處理器數量Np=128)增大到256×256×256(處理器數量Np=1024),測試結果見表7.結果表明,采用兩水平FS 預條件子的Newton–Krylov 求解器具有較好的并行性能,在線性迭代和計算時間方面均優于一水平FS 預條件子.

表6 算例一的強可擴展性(“–”表示未測試)

表7 算例一的弱可擴展性(“–”表示未測試)

同樣,研究算例二的一水平和兩水平FS 預條件子的強、弱可擴展性.表8為強可擴展性測試結果.測試中使用256×256×128 的網格進行模擬,處理器的數量從64 個增加到2048 個,比較非線性和線性迭代的數量以及相對應處理器數量的計算時間.表9為弱擴展性方面的并行性能結果.模擬從使用256 個處理器在256×256×128 的網格上計算開始,以1372 個處理器在448×448×2246網格上計算結束.結果表明,兩水平的FS 預條件子在線性迭代和計算時間方面優于一水平FS 方法.隨著處理器數量的增加,算法表現出良好的強、弱擴展性,CPR 預條件子優于加性FS 和乘性FS 預條件子.

表8 算例二的強可擴展性(“**”表示算法不收斂,“–”表示未測試)

表9 算例二的弱可擴展性(“**”表示算法不收斂,“–”表示未測試)

5 結論

本文研究用于裂縫型多孔介質非定常流動問題的場分裂預條件子.在區域分解框架下考慮幾種不同類型的場分裂預條件子,包括加性FS 預條件子、乘性FS 預條件子、舒爾補FS 預條件子和約束壓力殘差(CPR)預條件子.其中相應的子系統采用限制加性Schwarz(RAS)算法進行近似求解,RAS 的每個子塊通過LU 或ILU 分解得到.為了提高一水平FS 預條件子的計算效率和并行性能,在區域分解和多重網格方法的基礎上引入粗網格,設計兩水平FS 預條件子.最后在超級計算機上計算幾個基準案例,得到具有良好魯棒性和并行可擴展性的精確數值解.在對流動問題進行實驗后,發現兩水平的CPR 預條件子在線性迭代和計算時間方面對所研究的問題特別有效.本文所提出的多水平場分裂預條件子可以推廣到更復雜的流動問題.