多相流數(shù)值模擬中復(fù)雜地質(zhì)體網(wǎng)格剖分實(shí)現(xiàn)技術(shù)

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

（吉林大學(xué) 地下水資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)春130021）

0 引 言

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

1 PEBI網(wǎng)格

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

圖1 規(guī)則網(wǎng)格與PEBI網(wǎng)格Fig.1 Regular grid and PEBI grid

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

綜合目前的PEBI網(wǎng)格生成算法，其自動(dòng)化程度較高，沒有或較少有人工參與方面的生成方法，不利于研究者對(duì)具有多約束復(fù)雜邊界地質(zhì)體的網(wǎng)格剖分?；谶@種不足，本文利用可視化界面的交互式特點(diǎn)，提出運(yùn)用交互式布點(diǎn)的方法來生成復(fù)雜地質(zhì)條件下的PEBI網(wǎng)格，其算法簡(jiǎn)單，操作靈活、方便，可解決多相流模擬中復(fù)雜地質(zhì)體的空間網(wǎng)格離散問題。自動(dòng)化PEBI網(wǎng)格剖分技術(shù)的關(guān)鍵是剖分節(jié)點(diǎn)的生成，本文算法首先對(duì)研究區(qū)進(jìn)行合理布點(diǎn)，并進(jìn)行Delaunay［11］三角剖分，追蹤其對(duì)應(yīng)的泰森多邊形［12］，后分配高程，生成具有拓?fù)浣Y(jié)構(gòu)的PEBI網(wǎng)格。

圖2 PEBI網(wǎng)格模型Fig.2 Model of PEBI grid

2 網(wǎng)格生成中的約束條件

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

圖3 兩種情況下的PEBI網(wǎng)格Fig.3 PEBI grid in both cases

3 布點(diǎn)方法

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

3.1 布點(diǎn)原則

模擬區(qū)布置的點(diǎn)將作為生成PEBI網(wǎng)格的單元點(diǎn)（見圖4），點(diǎn)布置得合理與否將直接影響數(shù)學(xué)模型的求解計(jì)算規(guī)模、運(yùn)算時(shí)間與解的精度［13］。所以，在進(jìn)行網(wǎng)格剖分時(shí)需要注意以下幾方面問題：

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

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

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

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

3.2 布點(diǎn)方法

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

4 生成PEBI網(wǎng)格算法

布點(diǎn)完成后，接下來需要進(jìn)行三方面的工作：Delaunay三角剖分；查找邊界鈍角三角形，調(diào)整點(diǎn)布局；生成PEBI，以完成高質(zhì)量PEBI網(wǎng)格生成。

4.1 Delaunay三角剖分

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

4.2 邊界鈍角三角形

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

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

若余弦值等于零，則為直角三角形，小于零則為鈍角三角形?；诖嗽恚梢院苋菀渍业竭吔玮g角三角形（見圖5（b））。

圖5 PEBI網(wǎng)格剖分與修正Fig.5 PEBI grid subdivision and revision

4.3 PEBI網(wǎng)格生成

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

5 剖分實(shí)例與應(yīng)用

5.1 剖分實(shí)例

為了驗(yàn)證本文網(wǎng)格剖分方法的科學(xué)實(shí)用性，將算法耦合到作者前期開發(fā)的可視化建模軟件TOUGHVISUAL［18］平臺(tái)上進(jìn)行測(cè)試。

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

圖6 PEBI網(wǎng)格生成流程圖Fig.6 Flow chart of PEBI grid generation

圖7 多點(diǎn)／線約束PEBI網(wǎng)格生成Fig.7 PEBI grid generation under multi－point／line constrains

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

圖8（a）是一個(gè)點(diǎn)約束、一個(gè)線約束、二個(gè)區(qū)約束情況下的網(wǎng)格劃分結(jié)果。點(diǎn)約束區(qū)的網(wǎng)格剖分，從點(diǎn)約束處向外網(wǎng)格單元面積逐漸增大，以減緩數(shù)值計(jì)算時(shí)點(diǎn)約束區(qū)各物理場(chǎng)的變量梯度急劇變化的問題；線約束處強(qiáng)制要求單元不穿過線，并且單元密度從線約束處向外逐漸增大，以適應(yīng)線約束區(qū)向外的突變情況；區(qū)約束采用了兩種方式的網(wǎng)格單元類型，一種為矩形，如圖8（a）中的右上部分，另一種為多邊形網(wǎng)格單元，且單元密度由密到稀分布。這與點(diǎn)的布置直接相關(guān)。圖8（b）是通過組建相應(yīng)的二維網(wǎng)格拓?fù)浣Y(jié)構(gòu)，分配高程后生成的三維PEBI網(wǎng)格結(jié)果。

圖8 點(diǎn)、線、區(qū)約束PEBI網(wǎng)格生成Fig.8 PEBI grid generation under point，line and area constraints

5.2 應(yīng)用

二氧化碳地質(zhì)儲(chǔ)存是當(dāng)今研究的熱點(diǎn)問題之一。油氣藏是很好的二氧化碳地質(zhì)儲(chǔ)存場(chǎng)地，同時(shí)注入的二氧化碳可提高石油采收率［19］。本文以油田某一區(qū)塊地質(zhì)條件網(wǎng)格剖分建模為實(shí)例，進(jìn)行了本文方法的應(yīng)用。研究區(qū)為一背斜構(gòu)造，在注入井約500m 處有一斷層。圖9為運(yùn)用本文方法進(jìn)行網(wǎng)格剖分的過程，圖10為模擬區(qū)最后的網(wǎng)格剖分結(jié)果，圖11 為模擬計(jì)算的壓力場(chǎng)分布圖。

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

圖9 研究區(qū)三角網(wǎng)剖分示意圖Fig.9 Triangle mesh subdivision in the study area

圖10 研究區(qū)PEBI網(wǎng)格生成Fig.10 PEBI grid generation in study area

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

6 結(jié)束語

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

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

［2］凌建軍，吳敬軒.網(wǎng)格方向性對(duì)水驅(qū)油油藏?cái)?shù)值模擬結(jié)果的影響［J］.江漢石油學(xué)院學(xué)報(bào)，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網(wǎng)格的生成算法［J］.西南石油學(xué)院學(xué)報(bào)2006，28（2）：32－36.Xiang Zu－ping，Zhang Lie－h(huán)ui，Chen Zhong－h(huán)ua，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］蔡強(qiáng)，楊欽，孟憲海，等.二維PEBI網(wǎng)格的生成［J］.工程圖學(xué)學(xué)報(bào)，2005，26（2）：169－172.Cai Qiang，Yang Qin，Meng Xian－h(huán)ai，et al.Research on 2D PEBI grid generation［J］.Journal of Engineering Graphics，2005，26（2）：169－172.

［7］查文舒，李道倫，盧德唐，等.井間干擾條件下PEBI網(wǎng)格劃分研究［J］.石油學(xué)報(bào)，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網(wǎng)格的多分支水平井井型優(yōu)化研究［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網(wǎng)格精細(xì)油藏?cái)?shù)值模擬應(yīng)用研究［J］.石油學(xué)報(bào)，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］王代剛，侯健，邢學(xué)軍，等.基于前沿推進(jìn)的改進(jìn)型PEBI網(wǎng)格生成方法［J］.計(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］李海峰，吳冀川，劉建波，等.有限元網(wǎng)格剖分與網(wǎng)格質(zhì)量判定指標(biāo)［J］.中國(guó)機(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三角剖分與有限元網(wǎng)格自動(dòng)生成［J］.華中理工大學(xué)學(xué)報(bào)，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三角網(wǎng)的泰森多邊形生成算法研究［J］.長(zhǎng)江大學(xué)學(xué)報(bào)（自然科學(xué)版）理工卷，2007，4（1）：100－103.Liu Shao－h(huán)ua，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］李孟濤，單文文，劉先貴，等.超臨界二氧化碳混相驅(qū)油機(jī)理實(shí)驗(yàn)研究［J］.石油學(xué)報(bào)，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.