多相流數值模擬中復雜地質體網格剖分實現技術

楊艷林，靖 晶，楊志杰，許天福，王福剛

（吉林大學 地下水資源與環境教育部重點實驗室，長春130021）

0 引 言

多相流體廣泛存在于地下能源和資源開發以及廢物地質處置等多個領域［1］，如油藏開采、深部地熱能開發、核廢料地質處置、二氧化碳地質儲存等多個科學領域都涉及到氣相、液相、固相等多物理場耦合的科學問題。多相流數值模擬技術已經成為研究這些領域科學問題的重要工具，并得到了廣泛應用。在數值模擬中，對于復雜地質體的剖分建模技術一直以來都是影響數值模擬精準度的一個重要因素，不同的網格剖分方法對計算規模、計算結果產生很大影響。在多相流數值模擬領域，目前對于復雜地質體（如斷層、井周圍地質體）的剖分建模技術仍存在很大不足。本文基于前人的研究基礎，提出了基于布點法構建任意多邊形、任意約束的PEBI（Perpendicular bisection）多約束、交互式網格剖分實現技術與網格生成算法，有效解決了復雜地質體難以客觀刻畫的科學問題。

1 PEBI網格

在地質體網格剖分方法中，結構化網格作為最簡單的空間離散方法，應用廣泛，但在描述復雜地質條件時存在兩方面局限：①不能很好的反映復雜邊界情況，尤其是在處理井周圍和斷層等地質體時，結構化網格處理效果很不理想（見圖1）；凌建軍等［2］采用矩形網格來逼近任意形狀的邊界，但仍有較大的誤差；②存在不可忽略的網格取向效應，而且不利于網格加密。為了解決這一問題，有的學者［3－4］提出將PEBI網格應用到多相流模擬中。PEBI網格是一種非結構網格，利用了有限元劃分網格的靈活性，可以逼近任意油藏形狀，便于網格加密；同時，PEBI網格中相鄰網格塊的交界面垂直平分相應PEBI網格的連線，具有有限差分的剖分特點，最終得到的差分方程與笛卡爾差分方程相似，從而可更好地逼近流體流動形態，以獲得更加精確的數值解。

圖1 規則網格與PEBI網格Fig.1 Regular grid and PEBI grid

PEBI網格是一種局部正交網格，即任意兩個相鄰網格塊的交界面一定垂直平分相應PEBI網格節點的連線。國外學者在這方面做了大量的工作，Heinemann等［3］首次將PEBI網格應用到油藏模擬中，其后Palagi等［4］也將PEBI網格運用到實際油藏數值模型中，獲得了滿意的效果。國內在這方面的研究與應用雖起步較晚，但至今也取得了豐富的成果，向祖平等［5］先按點搜索的方法生成三角網，后再生成其對偶圖，連接其外接圓心生成PEBI網格（見圖2）；蔡強等［6］運用控制圓法生成了帶約束的PEBI網格；查文舒等［7］根據井間流線方程的特征，提出了井間干擾條件下的PEBI網格劃分算法。王星等［8］利用PEBI網格分別對多分支水平井井型優化進行研究；劉立明等［9］基于精細油藏數值模擬中出現的問題，提出了徑向網格和PEBI網格的混合PEBI算法；王代剛等［10］提出了基于前沿推進的改進型PEBI網格生成方法。

綜合目前的PEBI網格生成算法，其自動化程度較高，沒有或較少有人工參與方面的生成方法，不利于研究者對具有多約束復雜邊界地質體的網格剖分。基于這種不足，本文利用可視化界面的交互式特點，提出運用交互式布點的方法來生成復雜地質條件下的PEBI網格，其算法簡單，操作靈活、方便，可解決多相流模擬中復雜地質體的空間網格離散問題。自動化PEBI網格剖分技術的關鍵是剖分節點的生成，本文算法首先對研究區進行合理布點，并進行Delaunay［11］三角剖分，追蹤其對應的泰森多邊形［12］，后分配高程，生成具有拓撲結構的PEBI網格。

圖2 PEBI網格模型Fig.2 Model of PEBI grid

2 網格生成中的約束條件

在復雜的地質體網格剖分中，模擬區的非均質性、幾何形態（邊界斷層、地層尖滅等）的精細描述、巖石和流體物理性質的空間變化等都要求進行特殊的處理。本文將其綜合分成三類，概化為點、線、區等約束條件，如直井，可處理為點約束；水平井、斷層，可處理為線約束；其他一些具有區域性的約束則可處理為區約束。下文針對點、線、區約束，提供了一些布點方法，以達到事半功倍的效果。對于線約束，需要網格沿著折線分布，不允許出現跨過約束線的網格（見圖3）；同時網格的單元尺寸和形狀也需要做一定的控制，比如在井眼周圍，為了準確地反映井眼周圍流體流動特征，網格單元由小到大成放射狀分布。

圖3 兩種情況下的PEBI網格Fig.3 PEBI grid in both cases

3 布點方法

為了生成滿足幾何區域和數值條件要求的高質量網格，在網格點的布置過程中，節點的生成以及約束點的生成，需充分考慮給定幾何區域的形狀和大小、求解區域中解的變化情況和最終形成網格單元的質量，以獲得高質量高效率的數值模擬網格。

3.1 布點原則

模擬區布置的點將作為生成PEBI網格的單元點（見圖4），點布置得合理與否將直接影響數學模型的求解計算規模、運算時間與解的精度［13］。所以，在進行網格剖分時需要注意以下幾方面問題：

（1）確定合適的布點密度。在多相流（油藏）數值模擬中經常碰到單元網格應剖分得如何細致才能獲得合理結果的問題。當然，在模型運算前是很難回答這個問題的。一般的做法是先執行一個認為合理的網格密度進行試算，同時對于關鍵區域（如注入井）進行雙倍的網格加密，重新分析，并比較兩個模型結果。若結果幾乎相同，則網格剖分合理。若兩次結果相差顯著，則應繼續細化網格直到兩次獲得近似相等的結果。

（2）單元形狀與類型的選擇。剖分單元形狀可以是三角形、四邊形以及其他多邊形。剖分形狀盡量是正多邊形，這樣的計算精度更好些。同時也要考慮單元大小的過渡，如在加密區域，單元從稀到疏的逐漸過渡（見圖3）。

（3）網格方向（如井處的網格成放射狀）與流體的流動方向盡量保持一致。這樣的網格剖分更能反應流體流動特征，以提高計算精度。反之，則很容易在換算時，由于計算機的計算精度而引入誤差，誤差會在計算過程中傳遞下去，越來越大，使計算結果失真。

圖4 各種約束的布點方式Fig.4 Way of distribution point for different constraints

3.2 布點方法

綜合考慮實際布點的需要，本文主要討論了幾種常用的布點方法，如矩形布點法、環形布點法和推進布點法等。矩形布點主要是在給定的多邊形內生成矩形分布點，最后生成的PEBI單元也是矩形的。環形布點法是主要針對點約束而進行處理的一種方法，其生成的點環形向外擴散，以達到生成的PEBI網格方向與徑向流動方向相一致。推進布點法則主要用于處理多邊形布點，其在向內推進時，分布的點會少于前一次的點，反之，向外則增加（見圖4）。當然也可以先進行多邊形自動三角剖分，后進行PEBI網格劃分，這方面的參考文獻較多，如丁永祥［14］完成的自動三角剖分。

4 生成PEBI網格算法

布點完成后，接下來需要進行三方面的工作：Delaunay三角剖分；查找邊界鈍角三角形，調整點布局；生成PEBI，以完成高質量PEBI網格生成。

4.1 Delaunay三角剖分

Delaunay三角 剖 分［11］于1934 年 被 提 出，在數學、地理、工程等許多領域應用廣泛。目前常用的有逐點插入法、三角網生長法和分割合并法。插入點算法的步驟是：首先，定義一個包含所有節點的初始網格，最簡單的情形是單個三角形；向網格中插入一個節點，找出其外接圓包含此節點的所有三角形，刪除這些單元形成一個包含插入節點的空腔；將該插入節點與空腔的每條邊相連，形成新的三角形單元；上述的節點插入過程重復進行，直到全部節點插入完畢。

4.2 邊界鈍角三角形

網格剖分若出現鈍角三角形，將直接影響PEBI的網格質量，從而影響模型的計算精度。在邊界上的鈍角三角形會產生不在研究范圍內的PEBI網格，如圖5（a）所示，出現這種情況必須調整邊界上的點的位置，直至生成所有PEBI單元都在研究范圍內（見圖5（c））。邊界鈍角三角形的查找算法可根據余弦定理來進行判定。

由余弦定理可得角A（見圖5（d））的余弦值，見式（1）（其中a、b、c為邊的長度）：

若余弦值等于零，則為直角三角形，小于零則為鈍角三角形。基于此原理，可以很容易找到邊界鈍角三角形（見圖5（b））。

圖5 PEBI網格剖分與修正Fig.5 PEBI grid subdivision and revision

4.3 PEBI網格生成

泰森多邊形又叫馮洛諾伊圖［14］（Voronoi diagram），其應用廣泛，如北京的鳥巢設計就借鑒了泰森多邊形的思想，其原理是先求出每個三角形的外接圓圓心，后連接圓心，即可生成泰森多邊形，即PEBI網格。Atsuyuki等［15］、Brassel［16］、向祖 平等［5］、蔡 強 等［6］、劉 少 華 等［17］都 對 這 方 面 進行過研究。但是為了生成多相流（油藏）數值模擬可運行的網格以及三維空間離散，還需進行泰森多邊形的拓撲重構，其主要是基于查找共享三角形節點的PEBI單元來構建其拓撲結構，如圖2中角點A 周圍的O1、O2、O3、O4、O5、O6等組成的單元號。PEBI網格生成流程如圖6所示。

5 剖分實例與應用

5.1 剖分實例

為了驗證本文網格剖分方法的科學實用性，將算法耦合到作者前期開發的可視化建模軟件TOUGHVISUAL［18］平臺上進行測試。

圖7（a）為多井點約束（井網很密），流場存在相互干擾情況下PEBI網格剖分結果，點周圍為徑向網格，其他區域的網格為變尺度的網格。其特點是：距點約束（井）越遠，網格的尺度越大。區域中的3個約束點（丼）的徑向區域出現重疊，此時為了反映流動特征與真實的流動相一致，采用了干擾網格劃分，其他區域采用非干擾PEBI網格劃分。在布點時，應確定多點約束之間是否存在流場干擾，然后，在算法耦合的圖形界面上通過人機交互的控制方式進行網格剖分。

圖6 PEBI網格生成流程圖Fig.6 Flow chart of PEBI grid generation

圖7 多點／線約束PEBI網格生成Fig.7 PEBI grid generation under multi－point／line constrains

圖7 （b）為多線約束（如水平井與斷層交織等復雜情況）時采用多邊形與矩形混合的PEBI網格剖分結果。線約束區采用矩形網格，向外網格逐漸變稀，之后通過多邊形網格與約束區外矩形網格進行連接。約束區外剖分單元也可以是多邊形，在布點時，可以靈活控制。

圖8（a）是一個點約束、一個線約束、二個區約束情況下的網格劃分結果。點約束區的網格剖分，從點約束處向外網格單元面積逐漸增大，以減緩數值計算時點約束區各物理場的變量梯度急劇變化的問題；線約束處強制要求單元不穿過線，并且單元密度從線約束處向外逐漸增大，以適應線約束區向外的突變情況；區約束采用了兩種方式的網格單元類型，一種為矩形，如圖8（a）中的右上部分，另一種為多邊形網格單元，且單元密度由密到稀分布。這與點的布置直接相關。圖8（b）是通過組建相應的二維網格拓撲結構，分配高程后生成的三維PEBI網格結果。

圖8 點、線、區約束PEBI網格生成Fig.8 PEBI grid generation under point，line and area constraints

5.2 應用

二氧化碳地質儲存是當今研究的熱點問題之一。油氣藏是很好的二氧化碳地質儲存場地，同時注入的二氧化碳可提高石油采收率［19］。本文以油田某一區塊地質條件網格剖分建模為實例，進行了本文方法的應用。研究區為一背斜構造，在注入井約500m 處有一斷層。圖9為運用本文方法進行網格剖分的過程，圖10為模擬區最后的網格剖分結果，圖11 為模擬計算的壓力場分布圖。

由圖10可知，采用此種網格剖分方法，可以較好地反映地層和斷層的空間分布情況，證明了本文方法的有效性與正確性。由圖11可知，500 m 處的斷層對壓力場的重新分布產生了較大影響，符合實際情況。

圖9 研究區三角網剖分示意圖Fig.9 Triangle mesh subdivision in the study area

圖10 研究區PEBI網格生成Fig.10 PEBI grid generation in study area

圖11 壓力分布圖Fig.11 Pressure distribution

6 結束語

復雜地質體的剖分建模技術是影響數值模擬結果的重要環節，不同的網格處理方式，將直接影響計算穩定性、復雜程度和計算效率等。本文通過研究多相流數值模擬網格剖分的特點，提出復雜地質體網格剖分建模的技術和算法，并將該算法耦合到作者前期開發的多相流可視化建模界面軟件TOUGHVISUAL 上，為復雜地質體客觀、快速建模提供了技術實現，為深部能源和資源開發、核廢料與二氧化碳等廢物地質儲存中的多物理場多相流數值模擬提供了科學實用的工具。

［1］國麗萍，劉承婷，劉保君.石油工程多相流體力學［M］.北京：中國石化出版社，2011.

［2］凌建軍，吳敬軒.網格方向性對水驅油油藏數值模擬結果的影響［J］.江漢石油學院學報，1990，12（1）：40－45.Ling Jian－jun，Wu Jing－xuan.The influences of grid orientation on the results of reservoir simulation in oil reservoir drived by water［J］.Journal of Jiang－Han Petroleum Institute，1990，12（1）：40－45.

［3］Heinemann Z E，Brand C W.Gridding techniques in reservoir simulation［C］∥Proceedings of the First and Second International Forum on Reservoir Simulation，Alpbach，Austria，1989：339－426.

［4］Palagi C L，Aziz K.Use of voronoi grid in reservoir simulation［J］.SPE Advanced Technology Series，1994，2（2）：1－9.

［5］向祖平，張烈輝，陳中華，等.油藏任意約束平面域PEBI網格的生成算法［J］.西南石油學院學報2006，28（2）：32－36.Xiang Zu－ping，Zhang Lie－hui，Chen Zhong－hua，et al.Algorithm for constructing PEBI mesh of arbitrarily shaped and constrained planar domains in oil reservoir［J］.Journal of Southwest Petroleum Institute，2006，28（2）：32－36.

［6］蔡強，楊欽，孟憲海，等.二維PEBI網格的生成［J］.工程圖學學報，2005，26（2）：169－172.Cai Qiang，Yang Qin，Meng Xian－hai，et al.Research on 2D PEBI grid generation［J］.Journal of Engineering Graphics，2005，26（2）：169－172.

［7］查文舒，李道倫，盧德唐，等.井間干擾條件下PEBI網格劃分研究［J］.石油學報，2008，29（5）：742－746.Zha Wen－shu，Li Dao－lun，Lu De－tang，et al.PEBI grid division in inter－well interference area［J］.Acta Petrolet Sinca，2008，29（5）：742－746.

［8］王星，常鐵龍，馬艷，等.基于PEBI網格的多分支水平井井型優化研究［J］.鉆采工藝，2010，33（4）：45－48.Wang Xing，Chang Tie－long，Ma Yan，et al.Optimization of multi－lateral horizontal wells based on PEBI grid［J］.Drilling and Production Technology，2010，33（4）：45－48.

［9］劉立明，廖新維，陳欽雷.混合PEBI網格精細油藏數值模擬應用研究［J］.石油學報，2003，24（3）：64－67.Liu Li－ming，Liao Xin－wei，Chen Qin－lei.Usage of hybrid PEBI grid in fine reservoir simulation［J］.Acta Petrolei Sinca，2003，24（3）：64－67.

［10］王代剛，侯健，邢學軍，等.基于前沿推進的改進型PEBI網格生成方法［J］.計算物理，2012，29（5）：675－683.Wang Dai－gang，Hou Jian，Xing Xue－jun，et al.Modified PEBI grid generation method with advanceing front approach［J］.Chinese Journal of Computational Physics，2012，29（5）：675－683.

［11］Delaunay Boris.Sur la sphere vide［J］.Bulletin of the Academy of Sciences of the USSR，Classe des Sciences Mathematiques et Naturelles，1934（8）：793－800.

［12］Thiessen A H.Precipitation averages for large areas［J］.Monthly Weather Review，1911，39（7）：1082－1084.

［13］李海峰，吳冀川，劉建波，等.有限元網格剖分與網格質量判定指標［J］.中國機械工程，2012，23（3）：368－377.Li Hai－feng，Wu Ji－chuan，Liu Jian－bo，et al.Finite element mesh generation and decision criteria of mesh quality［J］.China Mechanical Engineering，2012，23（3）：368－377.

［14］丁永祥.約束Delaunay三角剖分與有限元網格自動生成［J］.華中理工大學學報，1995，23（6）：39－43.Ding Yong－xiang.The constrained delaunay triangulation and the automatic generation of finite element meshes［J］.Journal of Huazhong University of Science and Technology，1995，23（6）：39－43.

［15］Okabe A，Boots B，Sugihara K，et al.Spatial Tessellations：Concepts and Applications of Voronoi Diagrams［M］.2nd ed.Chichester：John Wiley ＆Sons，2008.

［16］Brassel K E，Reif D.A procedure to generate thiessen polygons［J］.Geographical Analysis，1979，11（3）：289－303.

［17］劉少華，羅小龍，何幼斌，等.基于Delauany三角網的泰森多邊形生成算法研究［J］.長江大學學報（自然科學版）理工卷，2007，4（1）：100－103.Liu Shao－hua，Luo Xiao－long，He You－bin，et al.Based on Delauany theissen polygon generation algorithm of triangulation［J］.Journal of Yangtze University（Natural Science Edition）Science and Engineering Volume，2007，4（1）：100－103.

［18］Yang Yan－lin，Xu Tian－fu，Wang Fu－gang，et al.A user friendly pre－processing and post－processing graphical interface graphical for toughreact transport of unsaturated groundwater and heat［C］∥TOUGH Symposium，Berkely，USA，2012：814－822.

［19］李孟濤，單文文，劉先貴，等.超臨界二氧化碳混相驅油機理實驗研究［J］.石油學報，2006，27（3）：80－83.Li Meng－tao，Shan Wen－wen，Liu Xian－gui，et al.Laboratory study on miscible oil displacement mechanism of supercritical carbon dioxide［J］.Acta Petrologica Sinica，2006，27（3）：80－83.