楊念 楊海建,* 邵柏強

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

1 引言

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

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

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

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

2 裂縫型多孔介質(zhì)雙孔雙滲模型

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

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

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

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

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

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

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

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

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

3 Field-Split 預(yù)條件子

本節(jié)介紹求解非線性代數(shù)系統(tǒng)的Newton-Krylov 算法,并在Newton-Krylov 算法中,構(gòu)建幾種基于物理和區(qū)域分解方法的的場分裂(FS)預(yù)條件子.

3.1 Newton-Krylov 算法

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

假設(shè)初始向量為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 算法中,求解器最重要的部分是選擇合適的預(yù)條件形式:

3.2 場分裂預(yù)條件子

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

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

加性FS 預(yù)條件子和乘性FS 預(yù)條件子分別為:

舒爾補FS 預(yù)條件子為:

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

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

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

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

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

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

3.3 限制加性的Schwarz 預(yù)條件子

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

限制加性Schwarz 預(yù)條件子定義為:

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

4 數(shù)值實驗

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

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?初始壓力隨深度的變化而變化,設(shè)定在10～11Mpa 之間?注水井和生產(chǎn)井的井底壓力分別為12Mpa和8Mpa.圖(b)–(d)為不同時刻裂縫壓力分布圖.

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

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

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

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

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

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

4.2 算例二:階梯型滲透率

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

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

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

表4 算例二中不同F(xiàn)S 預(yù)條件子性能的比較

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

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

4.3 算例三:SPE10 基準(zhǔn)算例

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

圖5 算例三中SPE10 基質(zhì)滲透率對數(shù)分布圖和孔隙度分布圖

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

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

4.4 算例四:隨機滲透率

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

圖6 算例四中隨機滲透率對數(shù)分布圖

4.5 并行測試

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

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

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

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

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

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

5 結(jié)論

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