2三角剖分算法介紹

易讀文庫易讀文庫易讀文庫2三角剖分算法介紹2.1Delaunay三角剖分理論基礎(chǔ)2.1.1VORONOI圖1908年，俄國數(shù)學(xué)家M.G.Voronoi提出了Voronoi圖[1]，它是計(jì)算幾何學(xué)上非常重要的工具之一，被廣泛應(yīng)用于各個(gè)領(lǐng)域中。Voronoi圖是由平面區(qū)域中連接兩鄰點(diǎn)的線段的中垂線所形成的區(qū)域。它又叫Dirichlet圖或Thiessen多邊形[2]。Voronoi圖是自然界中的宏觀物體和微觀物體以空間距離相互作用所形成的一種網(wǎng)格結(jié)構(gòu)，應(yīng)用的范圍相當(dāng)廣泛，特別是在計(jì)算幾何學(xué)領(lǐng)域的應(yīng)用[3]。Voronoi圖是一種關(guān)于平面或空間區(qū)域劃分的基礎(chǔ)數(shù)據(jù)結(jié)構(gòu)。100多年來，它被廣泛應(yīng)用在與幾何信息相關(guān)的各個(gè)領(lǐng)域。隨著計(jì)算機(jī)科學(xué)技術(shù)的不斷普及和發(fā)展，Voronoi圖的應(yīng)用領(lǐng)域也在不斷地?cái)U(kuò)大。對(duì)于Voronoi圖的應(yīng)用，以90年代以來的應(yīng)用更為突出。 Voronoi圖的本質(zhì)屬性是由空間實(shí)體幾何唯一確定的，不是通過其他方法強(qiáng)加上去的。從Voronoi圖的數(shù)字幾何角度來看，它是針對(duì)平面n個(gè)離散點(diǎn)而言的，把平面分為若干區(qū)域，每一個(gè)區(qū)域包含一個(gè)點(diǎn)，該點(diǎn)所在的區(qū)域就是到該點(diǎn)距離最近的點(diǎn)的集合。 設(shè)p1，p2是平面上兩點(diǎn)，L是線段p1p2的中垂線，L將平面分成兩部分LL和LR，位于LL內(nèi)的點(diǎn)pl具有特性：d(pl,p1)<d(pl,p2)，其中d(pl,pi)表示pl與pi之間的歐幾里得距離，i=1，2。這意味著，位于LL內(nèi)的點(diǎn)比平面上其他點(diǎn)更接近點(diǎn)p1，換句話說，LL內(nèi)的點(diǎn)是比平面上其他店更接近于p1的點(diǎn)的軌跡，記為V(p1)，如圖2.1所示。如果用H(p1,p2)表示半平面LL，而LR=H(p2,p1)，則有V(p1)=H(p1,p2)，V(p2)=H(p2,p1)。 給定平面上n個(gè)點(diǎn)的點(diǎn)集S，S={p1,p2,…,pn}。定義： V(pi)=即V(pi)表示比其他點(diǎn)更接近pi的點(diǎn)的軌跡是n-1個(gè)半平面的交，它是一個(gè)不多于n-1條邊的凸多邊形區(qū)域，稱為關(guān)于pi的Voronoi多邊形或關(guān)于pi的Voronoi域。對(duì)于S中的每個(gè)點(diǎn)都可以作一個(gè)Voronoi多邊形，這樣的n個(gè)Voronoi多邊形組成的圖稱為Voronoi圖。Voronoi多邊形的每條邊是S中某兩點(diǎn)的連線的垂直平分線，所有這樣的兩點(diǎn)連線構(gòu)成一個(gè)圖，稱為Voronoi圖的直線對(duì)偶圖。如圖2.2所示，虛線表示voronoi圖，實(shí)線表示其對(duì)偶圖。圖2.1平面中兩離散點(diǎn)圖2.2Voronoi圖 Voronoi圖有以下一些性質(zhì)：性質(zhì)1.n個(gè)點(diǎn)的點(diǎn)集S的Voronoi圖至多有2n-5個(gè)頂點(diǎn)和3n-6條邊。性質(zhì)2.每個(gè)Voronoi點(diǎn)點(diǎn)好是三條Voronoi邊的交點(diǎn)。性質(zhì)3.Voronoi圖的對(duì)偶圖實(shí)際上時(shí)點(diǎn)集的一種那個(gè)三角剖分，該三角剖分就是Delaunay三角剖分[4]。 性質(zhì)3實(shí)際上告訴我們可以用構(gòu)造Voronoi圖的方法來求平面點(diǎn)集的三角剖分，但實(shí)際上這種方法較少使用，因?yàn)樗惴ǖ男什缓谩?.1.2Delaunay三角剖分的基本概念1、基本概念（1）域分割 給定平面上的n個(gè)不相重的散亂數(shù)據(jù)點(diǎn)，對(duì)每個(gè)散亂數(shù)據(jù)點(diǎn)構(gòu)造一個(gè)域，使該域內(nèi)的任一點(diǎn)離此散亂點(diǎn)比離其他散亂點(diǎn)更近，這種域分割就是上節(jié)介紹的Voronoi圖，也稱Dirichlet域分割，由上節(jié)可知，域邊界其實(shí)就是連接兩個(gè)相鄰散亂點(diǎn)的直線的垂直平分線。 （2）Delaunay三角化 對(duì)平面上的散亂數(shù)據(jù)點(diǎn)進(jìn)行域分割后，將具有公共域邊界的散亂點(diǎn)對(duì)相連形成的三角剖分稱為Delaunay三角剖分。 （3）優(yōu)化 對(duì)三角網(wǎng)格進(jìn)行優(yōu)化，就是要使三角網(wǎng)格整體上盡量均勻，避免出現(xiàn)狹長三角形，也就是獲得Delaunay三角化。2、三角剖分優(yōu)化準(zhǔn)則 在三角剖分過程中，人們往往先用一種簡單的方法構(gòu)造散亂點(diǎn)的初始三角剖分，然后對(duì)其進(jìn)行優(yōu)化以獲得Delaunay剖分。優(yōu)化的方法取決于所采用的優(yōu)化準(zhǔn)則，平面三角剖分最常用的優(yōu)化準(zhǔn)則有Thiessen區(qū)域準(zhǔn)則、最小內(nèi)角最大準(zhǔn)則、圓準(zhǔn)則。Sibson證明了這三個(gè)準(zhǔn)則的等價(jià)性，并指出符合這三個(gè)準(zhǔn)則的三角剖分只有一個(gè)，即Delaunay三角剖分。最小內(nèi)角最大準(zhǔn)則對(duì)一個(gè)嚴(yán)格凸的四邊形進(jìn)行三角化時(shí)，有兩種選擇，最小內(nèi)角最大準(zhǔn)則就是要保證對(duì)角線兩側(cè)兩個(gè)三角形中的最小內(nèi)角為最大，如圖2.3所示。圖2.3對(duì)角線優(yōu)化圓準(zhǔn)則嚴(yán)格凸四邊形中的三個(gè)頂點(diǎn)確定一個(gè)圓，如果第四個(gè)頂點(diǎn)落在圓內(nèi)，則將第四個(gè)頂?shù)着c其相對(duì)的頂點(diǎn)相連，否則將另外兩個(gè)頂點(diǎn)相連，這個(gè)準(zhǔn)則稱為圓準(zhǔn)則。也就是說，符合圓準(zhǔn)則的三角剖分中，任一三角形的外接圓內(nèi)不應(yīng)該包含其他點(diǎn)，如圖2.4所示。圖2.4空外接圓準(zhǔn)則局部優(yōu)化局部優(yōu)化是指對(duì)任意一個(gè)凸四邊形的對(duì)角線，依據(jù)某種優(yōu)化準(zhǔn)則做交還測試后所得到的三角剖分。全局優(yōu)化當(dāng)三角形剖分T中每一條內(nèi)邊上的兩個(gè)三角形所形成的凸四邊形都滿足局部優(yōu)化標(biāo)準(zhǔn)時(shí)，稱該三角剖分T滿足全局優(yōu)化。退化前面已經(jīng)講過Delaunay三角剖分滿足圓準(zhǔn)則，即任一三角形的外接圓內(nèi)不能包含其他的點(diǎn)，如果規(guī)定三角剖分中任一三角形的外接圓內(nèi)和外接圓上都不能有其余頂點(diǎn)，則稱為標(biāo)準(zhǔn)的Delaunay三角剖分，在實(shí)際應(yīng)用中，往往存在四點(diǎn)共圓的情況。我們把四點(diǎn)或四點(diǎn)以上的散亂數(shù)據(jù)點(diǎn)共圓的情況稱為退化，相應(yīng)的三角剖分稱為準(zhǔn)Delaunay三角剖分。3、二維任意區(qū)域內(nèi)點(diǎn)集的Delaunay三角劃分概念及基本定理 將所有共邊的Voronoi多邊形的中心連接起來所形成的三角網(wǎng)稱為Delaunay三角剖分[5]。 在區(qū)域D內(nèi)的點(diǎn)集V的三角劃分（記為T(D,V)）是指將V中的點(diǎn)用互不相交的直線段連接起來，使得域D內(nèi)的每一個(gè)區(qū)域都是三角形，同時(shí)這些連線和三角形都在區(qū)域D內(nèi)，包含在區(qū)域D內(nèi)的邊稱為內(nèi)邊。 二維任意區(qū)域點(diǎn)集的Delaunay三角劃分[5]（簡記為DATA）是指所有的內(nèi)邊都滿足局部優(yōu)化的T(D,V)。4、基本定理 下面的定理是閔衛(wèi)東在文獻(xiàn)[6]中提出來的，其證明參見文獻(xiàn)[6]。定理2.1：設(shè)B，P點(diǎn)在AC邊的同側(cè)，如果角APC小于角ABC，則三角形ABC的外接圓內(nèi)不包含點(diǎn)P，參見圖2.5。圖2.5三角形外接圓圖2.6局部優(yōu)化邊AB定理2.2：設(shè)AB邊為三角形ABC和三角形ABE所共有，如果四邊形ACBE的內(nèi)角CBE≥180°或CAE≧180°，則AB邊是局部優(yōu)化的，見圖2.6。定理2.3：設(shè)三角形ABC是DATA（D，V）的一個(gè)三角形，P是V中的一點(diǎn)且不是ABC的三個(gè)頂點(diǎn)，如果線段PC在域D內(nèi)且它與AB邊有交點(diǎn)，則ABC的外接圓內(nèi)部不包含點(diǎn)P。定理2.4：（DATA的Circle準(zhǔn)則）三角剖分T（D，V）是DATA（D，V）的充要條件是，對(duì)于T（D，V）的任何一個(gè)三角形ABC，V中的任何點(diǎn)P，如果它與A，B，C相連的線段PA，PB，PC都在域D內(nèi)，則P不包含在三角形ABC的外接圓內(nèi)部。定理2.5：任何一個(gè)T（D，V）都可以通過有限次的局部優(yōu)化操作將其轉(zhuǎn)化為DATA（D，V）。定理2.6：DATA（D，V）在所有的T（D，V）中平均形態(tài)比最大。 以上幾個(gè)定理實(shí)際上作為Delaunay三角劃分呢算法的理論依據(jù)在許多以Delaunay原理為基礎(chǔ)的網(wǎng)格剖分算法中經(jīng)常用到。5、三角剖分時(shí)需要滿足的幾個(gè)約束條件 在進(jìn)行網(wǎng)格剖分過程中，生成三角形時(shí)，必須滿足一下三個(gè)約束條件：三角形相互之間是不相交的，即兩個(gè)三角形除端點(diǎn)外不應(yīng)該有別的交點(diǎn)，此條件稱為三角形相交約束。三角形相互之間是互不包含的，即任意一個(gè)三角形不能完全包含其他的三角形，此條件稱為三角形包含約束。三角形全部落在區(qū)域之內(nèi)，此條件稱為三角形合法約束。2.2Delaunay三角剖分算法介紹Delaunay三角剖分算法是最常用的一種平面點(diǎn)集的剖分算法，實(shí)際上，Delaunay三角剖分算法是基于Delaunay原理的一類三角剖分算法的總稱?？傮w來說，Delaunay三角剖分算法可分為以下幾種。2.2.1翻邊算法1977年，Lawson基于最小內(nèi)角最大優(yōu)化準(zhǔn)則提出了一種二維點(diǎn)集Delaunay三角剖分翻邊算法[7]。該算法的思想是：給定一個(gè)二維離散點(diǎn)集，首先對(duì)其進(jìn)行一次初始三角剖分，然后根據(jù)最小內(nèi)角最大化準(zhǔn)則判斷初始三角剖分中形成凸四邊形的共邊三角形是否滿足優(yōu)化準(zhǔn)則；如果不滿足，就交換四邊形的兩條對(duì)角線。一次循環(huán)直到所偶的三角形都滿足最小角最大優(yōu)化準(zhǔn)則為止。通過計(jì)算可以得到，翻邊算法的時(shí)間復(fù)雜度為O(n2)。根據(jù)前面的理論分析知道，運(yùn)用翻邊算法撲粉出來的所有三角形邊都是Delaunay三角剖分的邊，那么最終得到的三角網(wǎng)格一定是Delaunay三角網(wǎng)格。 目前，二維的翻邊算法發(fā)展比較完善，利用翻邊法對(duì)二維離散點(diǎn)集進(jìn)行Delaunay三角剖分原理比較簡單、得到的網(wǎng)格質(zhì)量也較好。1989年，BarryJoe基于空外接球準(zhǔn)則提出了一種三維Delaunay三角剖分的翻面法[8]。翻面法的原理跟翻邊法的原理是一樣的，只不過在三維空間利用翻面法會(huì)比較復(fù)雜，有待進(jìn)一步研究。2.2.2逐點(diǎn)插入算法隨著Delaunay三角剖分的不斷發(fā)展，大量剖分算法不斷被提出來。1981年，一種逐點(diǎn)插入算法被Bowyer和Watson提出來[9][10]。該算法的基本思想是：首先構(gòu)造點(diǎn)集的一個(gè)初始三角剖分，然后逐一地在當(dāng)前三角剖分中插入一個(gè)新點(diǎn)；在插入過程中，要根據(jù)Delaunay三角剖分空外接圓住著呢進(jìn)行三角網(wǎng)格優(yōu)化，直到整個(gè)點(diǎn)集為空集為止。 該算法首先構(gòu)造一個(gè)極大三角形，使得要剖分的點(diǎn)集中所有點(diǎn)都落在該三角形里面；其次，逐一地插入各個(gè)新點(diǎn)，每插入一個(gè)點(diǎn)，必須搜索位于外接圓內(nèi)部的點(diǎn)的三角形；然后，將這些三角形從三角形隊(duì)列中刪除，可以形成多邊形空腔；最后，將插入點(diǎn)和多邊形空腔連接起來，構(gòu)成若干個(gè)以插入點(diǎn)為共同頂點(diǎn)的新三角形。 定位新插入點(diǎn)的位置是逐點(diǎn)插入算法的一個(gè)關(guān)鍵步驟，直接關(guān)系到剖分網(wǎng)格的好壞和剖分速度的快慢。由于包含插入點(diǎn)的三角形的外接圓位于插入點(diǎn)的附近區(qū)域，只要能夠準(zhǔn)確定位新插入點(diǎn)在當(dāng)前三角剖分中的位置，便可快速檢索到外接圓包圍它的三角形。基于定位插入點(diǎn)位置的重要性，現(xiàn)在大量改進(jìn)的逐點(diǎn)插入算法主要都體現(xiàn)在定位新插入點(diǎn)的位置上。 由于逐點(diǎn)插入算法比較簡單、易懂，可以將其推廣到三維或者更高維空間的三角剖分。唯一不足的就是該算法時(shí)間復(fù)雜度較高，很多專家學(xué)者在研究算法的時(shí)候也會(huì)將該因素重點(diǎn)考慮進(jìn)去，從而設(shè)計(jì)出簡單、高效的三角剖分算法。2.2.3分割合并算法Hoey和Shamos提出了一種用于生成Voronoi圖的分割合并算法，由于Voronoi圖與Delaunay圖互為對(duì)偶圖，進(jìn)而可以間接得到Delaunay三角網(wǎng)格[11]。后來Dwyer、Katajaine、L.J.Guibas等人提出了用于直接生成Delaunay三角網(wǎng)格的分割合并算法。 分割合并算法主要根據(jù)數(shù)學(xué)遞歸思想，遞歸將點(diǎn)集分割成規(guī)模相當(dāng)?shù)膬刹糠肿蛹?，然后?duì)分割出來的點(diǎn)集進(jìn)行Delaunay三角剖分，再遞歸地將兩個(gè)相鄰Delaunay三角剖分進(jìn)行合并組合，從而生成了整個(gè)點(diǎn)集的Delaunay三角網(wǎng)格。改進(jìn)后的分割合并算法時(shí)間復(fù)雜度為O(nlogn)，目前最好的分割合并算法是Dewall算法。該算法不僅可以應(yīng)用于二維空間，還可對(duì)三維或更高維空間進(jìn)行三角剖分，雖然分割合并算法實(shí)現(xiàn)的時(shí)間效率很高，但它的編程卻非常復(fù)雜。 以上三種算法都是比較典型的Delaunay三角剖分算法，應(yīng)用范圍比較廣泛。除了這幾種經(jīng)典算法外，其實(shí)還有很多種其他的Delaunay三角剖分算法。在選擇三角剖分算法的時(shí)候，必須從難易程度、效率高低和貼近實(shí)際與否等方面選擇，其中逐點(diǎn)插入算法是目前最為簡單、最為流行的一種Delaunay三角剖分算法，該算法思想簡單、易懂，可推廣到維數(shù)更高空間的三角剖分。2.2.4Bowyer-Watson算法Bowyer-Watson算法是一種先形成粗原始網(wǎng)格，再逐步插入新點(diǎn)進(jìn)行細(xì)化的算法，是一種逐點(diǎn)插入算法。假定已有初始剖分為T，向已有網(wǎng)格內(nèi)插入新點(diǎn)的Bowyer-Watson算法如下： Step1：插入一個(gè)新點(diǎn)P到現(xiàn)有的Delaunay三角剖分中。 Step2：尋找并刪除所有外接圓包含P點(diǎn)的三角元素，形成一個(gè)Delaunay空腔。 Step3：連接P點(diǎn)和Delaunay腔壁上所有各點(diǎn)，形成新的Delaunay三角剖分。 應(yīng)用Bowyer-Watson算法進(jìn)行網(wǎng)格剖分的整個(gè)過程簡述如下：Step1：家里一個(gè)覆蓋整個(gè)計(jì)算區(qū)域的粗原始網(wǎng)格。Step2：采用Bowyer-Watson算法將邊界點(diǎn)逐點(diǎn)插入原始網(wǎng)格。邊界點(diǎn)由人工事先給定，并假設(shè)邊界點(diǎn)的分布式合理的。Step3：刪除計(jì)算區(qū)域以外的三角形，并確保邊界面的正確三角剖分（進(jìn)行拓?fù)湎嗳菪蕴幚恚?。Step4：用Bowyer-Watson算法逐步向計(jì)算區(qū)域內(nèi)插入新點(diǎn)（對(duì)于給定點(diǎn)集，按一定順序插入即可，若內(nèi)點(diǎn)是自動(dòng)生成的，則需要按一定的策略生成內(nèi)點(diǎn)，直到所有內(nèi)點(diǎn)都被插入或者網(wǎng)格達(dá)到一定的撲粉要求為止。Step5：對(duì)所生成的網(wǎng)格進(jìn)行拓?fù)湎嗳菪詸z查、網(wǎng)格光順等工作。3礦井涌水量地質(zhì)模型的實(shí)現(xiàn)3.1初始粗網(wǎng)格的建立礦井涌水量地質(zhì)模型的建立，首先是要將煤層頂?shù)装宓那骐x散成三角形網(wǎng)格。本文中采用逐點(diǎn)插入法作為主要的算法。 逐點(diǎn)插入法的第一步是要建立一個(gè)包含所有計(jì)算域頂點(diǎn)的初始粗網(wǎng)格。常用的建立初始網(wǎng)格的方法是包容盒法，一般可建立兩種形式的包容盒，一種是建立矩形的包容盒[12]，另一種是三角形包容盒[13]。無論是給定點(diǎn)集或者是給定區(qū)域的三角剖分均可采用該方法。如果只給出邊界區(qū)域，則內(nèi)點(diǎn)需要自動(dòng)產(chǎn)生。本文所建立的三角剖分是給定區(qū)域邊界，自動(dòng)生成內(nèi)點(diǎn)。這種方法在使用過程中需要對(duì)生成的三角形進(jìn)行邊界拓?fù)湎嗳菪詸z測。因?yàn)橥負(fù)湎嗳菪詸z測的步驟，本文給出一種無需建立三角形粗網(wǎng)格或矩形粗網(wǎng)格的方法，直接將給定區(qū)域的多邊形劃分成三角形初始粗網(wǎng)格。 該方法首先將給定邊界按順時(shí)針或者逆時(shí)針的方向?qū)⑦吔珥旤c(diǎn)進(jìn)行排列，然后計(jì)算每一個(gè)頂點(diǎn)的角度大小，作為判斷頂點(diǎn)角度的指標(biāo)。隨后每次選擇最小的頂點(diǎn)，將其與相鄰的頂點(diǎn)連成一個(gè)三角形，重新計(jì)算該頂點(diǎn)兩側(cè)的頂點(diǎn)角度大小，然后再選擇最小角進(jìn)行處理。 由于多邊形并非是完全的凸多邊形，邊界上很可能會(huì)出現(xiàn)凹頂點(diǎn)。在這種情況下，需要判斷頂點(diǎn)的凹凸性。 本文使用矢量面積法[14]來判斷多邊形頂點(diǎn)的凹凸性。如果多邊形邊界以逆時(shí)針方向排列，那么為了判斷某一頂點(diǎn)的凹凸性，可以用與該頂點(diǎn)相鄰的兩個(gè)頂點(diǎn)，共三個(gè)頂點(diǎn)來判斷該頂點(diǎn)的凹凸性。計(jì)算這三個(gè)頂點(diǎn)形成的三角形的矢量面積，如果結(jié)果為正，那么該頂點(diǎn)為凸頂點(diǎn)，如果為負(fù)，那么該頂點(diǎn)為凹頂點(diǎn)。計(jì)算公式可簡單描述如下： 設(shè)三個(gè)頂點(diǎn)為：A、B、C，那么： 判斷完頂點(diǎn)的凹凸性之后，需要計(jì)算頂點(diǎn)的角度（弧度）?？梢允褂孟噜徣齻€(gè)頂點(diǎn)形成的兩個(gè)向量來計(jì)算頂點(diǎn)的角度。、的數(shù)量積可以計(jì)算出兩向量的夾角，再結(jié)合頂點(diǎn)凹凸性的判斷，加上或者減去180°便是該頂點(diǎn)的角度。Delaunay剖分的初始化步驟可用以下算法描述：算法3.1Step1：判斷頂點(diǎn)凹凸性，計(jì)算每一個(gè)頂點(diǎn)的角度。Step2：從頂點(diǎn)數(shù)組中找到角度最小的頂點(diǎn)，連接該頂點(diǎn)兩側(cè)的頂點(diǎn)，形成一個(gè)三角形。并計(jì)入三角形數(shù)組。Step3：重新計(jì)算該頂點(diǎn)兩側(cè)新頂點(diǎn)的角度，并從頂點(diǎn)數(shù)組中刪去該頂點(diǎn)。Step4：從新的頂點(diǎn)數(shù)組中找到新的最小角度頂點(diǎn)，重復(fù)step2至step3。直到頂點(diǎn)數(shù)組中的頂點(diǎn)個(gè)數(shù)等于3個(gè)。Step5：連接最后三個(gè)頂點(diǎn)形成一個(gè)三角形，計(jì)入三角形數(shù)組。按照以上算法，最后會(huì)形成一個(gè)初始的粗三角形網(wǎng)格，或者稱為多邊形的一個(gè)粗三角劃分。如圖3.1，3.2。圖3.1任意多邊形區(qū)域圖3.2初始網(wǎng)格劃分3.2簡單多邊形的Delaunay剖分在建立了多邊形的初始粗網(wǎng)格以后，便可以開始向計(jì)算區(qū)域內(nèi)插入點(diǎn)進(jìn)行Delaunay三角剖分。本文采用自動(dòng)布點(diǎn)技術(shù)向計(jì)算區(qū)域內(nèi)插入點(diǎn)。采用的方法是：選擇面積最大三角形的最長邊的中點(diǎn)最為插入點(diǎn)，建立三角網(wǎng)格。這樣需要設(shè)置一個(gè)面積的控制量，以便控制剖分的程度。即三角形的個(gè)數(shù)或者三角網(wǎng)格的疏密程度。這是一個(gè)可以改變的量，可以方便地控制三角網(wǎng)格的大小。當(dāng)插入一個(gè)點(diǎn)以后，需要搜索整個(gè)三角形數(shù)組，進(jìn)行空外接圓檢測，如果新插入點(diǎn)落在該三角形的外接圓內(nèi)（包含邊界），那么該三角形被添加到需要優(yōu)化的三角形數(shù)組中，所有三角形檢測完后，符合要求的三角形會(huì)形成一個(gè)多邊形空腔，刪除這些三角形的公共邊，連接插入點(diǎn)與多邊形空腔的各個(gè)頂點(diǎn)，形成新的多個(gè)三角形，加入三角形數(shù)組，刪除被合并的三角形。這樣一直到三角形面積符合控制要求結(jié)束。便完成了簡單多邊形的Delaunay剖分。算法可描述如下：算法3.2 Step1：設(shè)定三角形面積控制單元，計(jì)算三角形數(shù)組中所有三角形的面積大小。 Step2：選擇面積最大的三角形的最長邊作為新的插入點(diǎn)，并將該點(diǎn)加到頂點(diǎn)數(shù)組的后面。 Step3：搜索整個(gè)三角形數(shù)組，找出所有插入點(diǎn)落在三角形外接圓中的三角形，計(jì)入優(yōu)化三角形數(shù)組，并將三角形的邊計(jì)入邊數(shù)組，并且不能有重復(fù)邊，若有重復(fù)邊，便要從數(shù)組中刪除該邊。Step4：連接插入點(diǎn)與邊數(shù)組中每條邊的兩個(gè)頂點(diǎn)，形成新的三角形，加入三角形數(shù)組。形成新的三角形數(shù)組。 Step5：刪除優(yōu)化三角形數(shù)組中的三角形，更新整個(gè)三角形數(shù)組，完成一次點(diǎn)的插入。 Step6：重復(fù)step2、step3、step4，直到三角形的面積符合控制要求，結(jié)束點(diǎn)的插入，這樣就完成了簡單多邊形的Delaunay剖分。剖分實(shí)例如圖3.3。圖3.3Delaunay剖分在剖分的過程中需要注意幾個(gè)問題。如果插入點(diǎn)落在多邊形邊界上，便會(huì)形成同一直線上的兩條不同的邊被同時(shí)計(jì)入邊數(shù)組，當(dāng)連接三角形的時(shí)候，便會(huì)出現(xiàn)三個(gè)頂點(diǎn)落在一條直線上的三角形，所以需要對(duì)這種三角形進(jìn)行特別的處理。不然會(huì)出現(xiàn)錯(cuò)誤。處理方法很簡單，在連接三角形的時(shí)候，需要計(jì)算插入點(diǎn)與兩個(gè)頂點(diǎn)連接形成的直線的斜率，比較這兩條直線的斜率，如果斜率相等，那么表示這三個(gè)點(diǎn)處在同一條直線上，這樣便不能連接三角形。算法表示如下：算法3.3 Step1：計(jì)算邊數(shù)組中將要連接的那條邊的兩個(gè)頂點(diǎn)與插入點(diǎn)形成的直線的斜率。 Step2：判斷兩條直線的斜率。 Step3：如果斜率相等，那么那條邊的連接性置為0，如果斜率不等那么直線的連接性置為1。Step4：連接邊的兩個(gè)頂點(diǎn)與插入點(diǎn)，形成一個(gè)三角形，計(jì)入三角形數(shù)組。在計(jì)算斜率的時(shí)候，由于數(shù)據(jù)類型的限制，三角形的斜率不一定嚴(yán)格相等，所以當(dāng)三點(diǎn)一線的時(shí)候斜率會(huì)有微小的誤差，這一誤差需要用一常量來作為判斷斜率是否相等的限制，可以將這一誤差設(shè)定為1°，用來幫助判斷斜率。3.3用VB實(shí)現(xiàn)的Delaunay剖分算法 本文使用VisualBasic6.0作為開發(fā)環(huán)境來實(shí)現(xiàn)Delaunay剖分算法。分為四個(gè)模塊：公共數(shù)據(jù)模塊初始化模塊剖分模塊數(shù)據(jù)輸出模塊3.3.1公共數(shù)據(jù)模塊 公共數(shù)據(jù)模塊定義有程序中使用的數(shù)據(jù)類型，包括：頂點(diǎn)類，存儲(chǔ)有頂點(diǎn)的x、y、z坐標(biāo)數(shù)據(jù)，以及頂點(diǎn)是否是邊界格點(diǎn)的判斷變量。三角形類，存儲(chǔ)有每個(gè)三角形的三個(gè)頂點(diǎn)編號(hào)，編號(hào)既是三角形頂點(diǎn)在頂點(diǎn)數(shù)組中的序號(hào)。臨時(shí)頂點(diǎn)類，該數(shù)據(jù)類型主要用于初始化時(shí)存儲(chǔ)臨時(shí)頂點(diǎn)數(shù)據(jù)，除了包括x、y、z坐標(biāo)外，還有弧度值與該頂點(diǎn)在原頂點(diǎn)數(shù)組中的序號(hào)。邊類，該數(shù)據(jù)類主要用于剖分時(shí)存儲(chǔ)被優(yōu)化三角形的邊信息，數(shù)組中的邊將成為新生成的三角形的邊。模塊中還定義有全局頂點(diǎn)數(shù)組和全局三角形數(shù)組，全局頂點(diǎn)個(gè)數(shù)和全局三角形個(gè)數(shù)。全局頂點(diǎn)數(shù)組是一個(gè)只增不減的數(shù)組，頂點(diǎn)數(shù)會(huì)隨著剖分的進(jìn)行而不斷增加，且數(shù)據(jù)不會(huì)改變。全局三角形數(shù)組是一個(gè)動(dòng)態(tài)數(shù)組，數(shù)組元素的個(gè)數(shù)和元素的取值都會(huì)隨著程序的運(yùn)行而改變。初始化時(shí)已事先為兩個(gè)數(shù)組定義了元素個(gè)數(shù)?？呻S剖分的細(xì)分程度修改。3.3.2初始化模塊 初始化模塊用于對(duì)多邊形區(qū)域進(jìn)行初始的網(wǎng)格化。形成初始的三角形存入三角形數(shù)組。該模塊中先定義臨時(shí)頂點(diǎn)數(shù)組，該臨時(shí)頂點(diǎn)數(shù)組繼承全局頂點(diǎn)數(shù)組的全部數(shù)據(jù)。并存儲(chǔ)每一個(gè)頂點(diǎn)在全局頂點(diǎn)數(shù)組中的序號(hào)。 隨后計(jì)算每一個(gè)頂點(diǎn)的角度，存入臨時(shí)頂點(diǎn)數(shù)組，使用算法3.1判斷頂點(diǎn)的凹凸性并計(jì)算頂點(diǎn)的角度。經(jīng)測試矢量面積法能判斷頂點(diǎn)的凹凸性，由此可以計(jì)算每一個(gè)頂點(diǎn)的角度大小，在程序中，角度值轉(zhuǎn)換為弧度值計(jì)算。 該模塊中包括判斷頂點(diǎn)凹凸性的函數(shù)，初始化網(wǎng)格的過程和繪制網(wǎng)格的過程。這三個(gè)函數(shù)或過程在初始化過程中被調(diào)用。 初始化網(wǎng)格的過程中，先找到頂點(diǎn)角度最小的角，將該頂點(diǎn)兩側(cè)的頂點(diǎn)連接形成一個(gè)三角形，插入全局三角形數(shù)組。并在臨時(shí)頂點(diǎn)數(shù)組中刪除該頂點(diǎn)。重新計(jì)算該頂點(diǎn)兩側(cè)的頂點(diǎn)角度大小，并更新臨時(shí)頂點(diǎn)數(shù)組。在計(jì)算的過程中，需要考慮序號(hào)為1、2、n-1、n這四個(gè)頂點(diǎn)。因?yàn)閿?shù)組數(shù)據(jù)類型的限制，當(dāng)最小角頂點(diǎn)序號(hào)是這四個(gè)值時(shí)，頂點(diǎn)兩側(cè)的頂點(diǎn)序號(hào)并不連續(xù)，所以需要分開計(jì)算。用判斷語句分別判斷這四種情況，分別計(jì)算頂點(diǎn)兩側(cè)的頂點(diǎn)角度。其余情況可用循環(huán)語句計(jì)算。執(zhí)行完初始化網(wǎng)格過程，全局三角形數(shù)組包含最初的三角形網(wǎng)格，全局頂點(diǎn)數(shù)組不變。調(diào)用繪制網(wǎng)格過程繪制三角形。3.3.3剖分模塊該模塊用于將多邊形剖分成需要的Delaunay三角形網(wǎng)格。模塊首先建立邊數(shù)組，用于存儲(chǔ)需要連接的邊信息。首先計(jì)算現(xiàn)有三角形的面積，找到面積最大的三角形。然后取該三角形最長邊的中點(diǎn)作為新的插入點(diǎn)，將插入點(diǎn)計(jì)入全局頂點(diǎn)數(shù)組。搜索整個(gè)全局三角形數(shù)組，找到插入點(diǎn)落在其外接圓中的三角形，并記錄該三角形在全局三角形數(shù)組中的序號(hào)，將其三條邊送入邊數(shù)組等待優(yōu)化。在將邊送入邊數(shù)組的過程中，需要判斷邊是否已經(jīng)存在邊數(shù)組中，如果已經(jīng)存在則要在邊數(shù)組中刪除該邊，如果不存在則存入邊數(shù)組。搜索完全部三角形后，形成了一個(gè)記錄需要優(yōu)化三角形的序號(hào)數(shù)組和邊數(shù)組，邊數(shù)組中的邊形成一個(gè)連續(xù)的多邊形空腔，插入點(diǎn)可能位于空腔內(nèi)部，或者空腔的某一條邊上。如果位于邊上，則需要檢測插入點(diǎn)落在哪條邊上，而這條邊則不能作為新三角形的邊。使用插入點(diǎn)與邊兩個(gè)頂點(diǎn)的連線的斜率來判斷插入點(diǎn)是否落在邊上，這樣便可除去三點(diǎn)一線的三角形。按順序連接插入點(diǎn)與邊數(shù)組中的每條邊，形成新的數(shù)個(gè)三角形，將這些三角形送入全局三角形數(shù)組。如果符合要求的三角形只有一個(gè)，則要專門計(jì)算插入三角形的個(gè)數(shù)，形成兩個(gè)新的三角形，刪除一個(gè)舊的三角形。隨后檢測需要?jiǎng)h除的三角形，從數(shù)組尾端開始向前逐個(gè)刪除需要?jiǎng)h除的三角形。這樣便完成了一個(gè)點(diǎn)在多邊形區(qū)域中的插入，更新了全局頂點(diǎn)數(shù)組和全局三角形數(shù)組。重新計(jì)算數(shù)組中每個(gè)三角形的面積，找到面積最大的三角形，如果面積值大于某一設(shè)定值，則繼續(xù)向其中插入點(diǎn)，直到滿足要求為止。這樣便完成了多邊形區(qū)域的Delaunay剖分。3.3.4數(shù)據(jù)輸出模塊 該模塊主要用于頂點(diǎn)數(shù)據(jù)和三角形數(shù)據(jù)的輸出，以便在隨后的計(jì)算中使用。將需要的數(shù)據(jù)輸出到一個(gè)文本文件中，留待后用。程序界面如圖3.4所示：圖3.4程序運(yùn)行界面3.4對(duì)剖分算法的評(píng)價(jià)該算法可以將簡單的多邊形剖分成Delaunay三角網(wǎng)格。并且算法省去了建立初始大三角形或者矩形的步驟，直接將多邊形劃分成初始粗網(wǎng)格。然后進(jìn)行剖分。這樣就省去了進(jìn)行邊界拓?fù)湎嗳菪詸z查的步驟。和其他算法相比，該算法在時(shí)間復(fù)雜度上并無優(yōu)勢，剖分本質(zhì)上是逐點(diǎn)插入法，時(shí)間復(fù)雜度較高。 該算法的一個(gè)缺點(diǎn)是，當(dāng)三角形長邊位于邊界上時(shí)，不能保證該三角形是銳角三角形，也是就該三角形有可能是鈍角三角形。目前尚無有效的辦法解決此問題，這也是選擇布點(diǎn)策略時(shí)產(chǎn)生的結(jié)果。3.5網(wǎng)格節(jié)點(diǎn)編號(hào)優(yōu)化3.5.1節(jié)點(diǎn)編號(hào)優(yōu)化算法有限差分分析時(shí)經(jīng)常會(huì)碰到具有稀疏矩陣剛度的方程式。減少稀疏矩陣的帶寬可以減少有限差分法的計(jì)算時(shí)間，和減少占用計(jì)算機(jī)的存儲(chǔ)空間。在有限差分法的計(jì)算過程中，可通過節(jié)點(diǎn)重新編號(hào)來減少帶寬。對(duì)于節(jié)點(diǎn)重編號(hào)的研究較多，最典型的有：Cuthill和Mckee[15]基于圖論的節(jié)點(diǎn)分層編號(hào)方法，可在一定程度上減少剛度矩陣的帶寬，這個(gè)方法后來被稱為CM方法。Gibbs，Poole和Stockmeyer[16]提出了更加有效的基于圖論的減少稀疏矩陣帶寬和外形的方法，此方法后來被稱為GPS方法。Akhras和Dhatt[17]通過計(jì)算節(jié)點(diǎn)商方法來縮小帶寬，該方法編程方便，此方法后來被稱為AD方法。歐陽興[18]提出了前沿法與矩形法。 CM方法和GPS方法，采用了復(fù)雜的圖論方法[19]，沒有估算算法本身的復(fù)雜度。實(shí)現(xiàn)比較困難，如果考慮節(jié)點(diǎn)編號(hào)算法本身時(shí)間消耗大，就會(huì)降低節(jié)點(diǎn)編號(hào)算法的效率。AD法雖然實(shí)現(xiàn)比較方便，但是在網(wǎng)格重劃和網(wǎng)格加密時(shí)節(jié)點(diǎn)編號(hào)優(yōu)化速度慢。前沿法和矩形法對(duì)于板料單元和平面二維單元特別適用。 黃志超，包忠詡，周天瑞[20]等從分析最佳最簡單的矩形出發(fā)，總結(jié)出一種新的節(jié)點(diǎn)重編號(hào)方法，是一種改進(jìn)的AD法。 如圖3.5是最典型的最佳二維網(wǎng)格節(jié)點(diǎn)編號(hào)，由此來分析其特性，總結(jié)出一般規(guī)律。圖3.5二維網(wǎng)格15節(jié)點(diǎn)優(yōu)化 該15節(jié)點(diǎn)序號(hào)是最佳的編號(hào)法，它的最優(yōu)同單元節(jié)點(diǎn)編號(hào)的最大差值為4。由該圖的編號(hào)法可總結(jié)出以下規(guī)則：節(jié)點(diǎn)商是從小到大升序排列的；最小最大節(jié)點(diǎn)編號(hào)和也是從小到大升序排列的；最大節(jié)點(diǎn)基本上是按升序排列的，不過有編號(hào)重復(fù)的情況。 由此可總結(jié)出一種新的節(jié)點(diǎn)重編號(hào)方法。算法3.4：Step1： 計(jì)算出每個(gè)節(jié)點(diǎn)的相關(guān)單元節(jié)點(diǎn)編號(hào)總和；Step2： 計(jì)算出該節(jié)點(diǎn)所相關(guān)單元的最大節(jié)點(diǎn)編號(hào)與最小節(jié)點(diǎn)編號(hào)；Step3： 計(jì)算出最小節(jié)點(diǎn)編號(hào)與最大節(jié)點(diǎn)編號(hào)之和；Step4：計(jì)算出該節(jié)點(diǎn)的節(jié)點(diǎn)商，將該節(jié)點(diǎn)的相關(guān)單元節(jié)點(diǎn)編號(hào)總和除以該節(jié)點(diǎn)的相關(guān)單元個(gè)數(shù)所得的值為節(jié)點(diǎn)商；Step5： 把節(jié)點(diǎn)商作為該節(jié)點(diǎn)重編號(hào)的依據(jù)，節(jié)點(diǎn)商小的編號(hào)也小，節(jié)點(diǎn)商大的編號(hào)也大，如果節(jié)點(diǎn)商相等則把最小最大節(jié)點(diǎn)編號(hào)和作為重編號(hào)的依據(jù)；step6： 重復(fù)（1）到（5）直到帶寬不能再減小為止。經(jīng)計(jì)算該算法能減少50%左右的帶寬。圖3.6為一495個(gè)節(jié)點(diǎn)的網(wǎng)格形成的稀疏矩陣每一行的帶寬分布圖。圖3.6節(jié)點(diǎn)帶寬分布圖3.5.2節(jié)點(diǎn)編號(hào)優(yōu)化算法的實(shí)現(xiàn) 本文使用VB語言實(shí)現(xiàn)節(jié)點(diǎn)編號(hào)算法。輸入數(shù)據(jù)為剖分完成后的節(jié)點(diǎn)數(shù)組和三角形數(shù)組。該算法對(duì)應(yīng)程序的計(jì)算水頭模塊。 在公共數(shù)據(jù)模塊中定義一新的數(shù)據(jù)類OptVertex。該數(shù)據(jù)類包含節(jié)點(diǎn)在頂點(diǎn)數(shù)組中的序號(hào)number，節(jié)點(diǎn)商N(yùn)odeTrade和最大最小節(jié)點(diǎn)編號(hào)之和maxmin三個(gè)數(shù)據(jù)。在計(jì)算水頭模塊定義一該類型的數(shù)組作為節(jié)點(diǎn)優(yōu)化數(shù)組，命名為OptVer。數(shù)組長度和全局頂點(diǎn)數(shù)組相同。 算法首先將全局頂點(diǎn)數(shù)組的序號(hào)存入節(jié)點(diǎn)優(yōu)化數(shù)組，然后開始循環(huán)。計(jì)算節(jié)點(diǎn)優(yōu)化數(shù)組的每一個(gè)節(jié)點(diǎn)，搜索全局三角形數(shù)組如果三角形中包含該節(jié)點(diǎn)的編號(hào)，便將該節(jié)點(diǎn)相關(guān)單元數(shù)加1，并根據(jù)三角形的三個(gè)頂點(diǎn)在全局頂點(diǎn)數(shù)組中的序號(hào)從節(jié)點(diǎn)優(yōu)化數(shù)組中找到這三個(gè)頂點(diǎn)，并將其在節(jié)點(diǎn)優(yōu)化數(shù)組中對(duì)應(yīng)的序號(hào)累加存入一變量中。同時(shí)要記錄相關(guān)單元中節(jié)點(diǎn)序號(hào)的最大值與最小值。三角形數(shù)組搜索完成后，便得到了該節(jié)點(diǎn)的相關(guān)單元數(shù)、相關(guān)單元節(jié)點(diǎn)編號(hào)總和、最大節(jié)點(diǎn)編號(hào)和最小節(jié)點(diǎn)編號(hào)。接著計(jì)算該節(jié)點(diǎn)的節(jié)點(diǎn)商和最大最小節(jié)點(diǎn)編號(hào)之和，分別存入該節(jié)點(diǎn)的對(duì)應(yīng)變量中。這樣依次計(jì)算所有節(jié)點(diǎn)優(yōu)化數(shù)組中所有節(jié)點(diǎn)的相關(guān)數(shù)據(jù)存入數(shù)組中。 對(duì)計(jì)算完成的節(jié)點(diǎn)進(jìn)行排序，使用選擇排序法，將第一個(gè)位置的節(jié)點(diǎn)和后面所有的節(jié)點(diǎn)做比較，節(jié)點(diǎn)商大的序號(hào)也大，節(jié)點(diǎn)商小的序號(hào)也小。如果節(jié)點(diǎn)商相等，比較節(jié)點(diǎn)的最大最小節(jié)點(diǎn)編號(hào)之和，小的序號(hào)也小，大的序號(hào)也大。然后比較第二個(gè)位置的節(jié)點(diǎn)和后面的所有節(jié)點(diǎn)，依次排序下去。這樣就完成了一次節(jié)點(diǎn)編號(hào)的優(yōu)化。 接著計(jì)算每一個(gè)節(jié)點(diǎn)的帶寬，初始帶寬為節(jié)點(diǎn)的個(gè)數(shù)。然后計(jì)算每一個(gè)節(jié)點(diǎn)在稀疏矩陣中的帶寬，計(jì)算方法為，從三角形數(shù)組中找到包含該節(jié)點(diǎn)的三角形，計(jì)算每一個(gè)節(jié)點(diǎn)和該節(jié)點(diǎn)之差的絕對(duì)值，也就是在稀疏矩陣的一行中包含數(shù)據(jù)的元素到該行對(duì)角線的距離的絕對(duì)值。計(jì)算出最大的帶寬，然后和初始帶寬比較，如果比初始帶寬小，說明稀疏矩陣帶寬減小了，應(yīng)當(dāng)進(jìn)行下一次節(jié)點(diǎn)編號(hào)優(yōu)化。將當(dāng)前帶寬作為下一次比較的帶寬，繼續(xù)循環(huán)重復(fù)以上步驟對(duì)節(jié)點(diǎn)編號(hào)進(jìn)行優(yōu)化排序。直到當(dāng)前帶寬不再改變?yōu)橹埂?這樣就完成了節(jié)點(diǎn)編號(hào)的優(yōu)化。3.6三角形網(wǎng)格有限差分法3.6.1三角形均衡網(wǎng)格差分方程的建立三角形網(wǎng)格差分方程是計(jì)算網(wǎng)格中格點(diǎn)水頭的重要步驟。也是有限差分方法的核心。差分方程形成的系數(shù)矩陣是主要的數(shù)據(jù)結(jié)構(gòu)。 考慮以格點(diǎn)i為中心的多邊形網(wǎng)格Di作為均衡區(qū)（如圖3.7），下面可根據(jù)達(dá)西定律和水均衡原理建立格點(diǎn)i的差分方程[21]。圖3.7多邊形網(wǎng)格均衡區(qū) 首先求出單位時(shí)間內(nèi)通過的周邊流入內(nèi)的水量。 由圖3.7可見，我們可以先求出通過三角形ijk內(nèi)的邊和流入內(nèi)的水量。由達(dá)西定律有：（3.1）式中：、為流段ij和ik的平均導(dǎo)水系數(shù)；、、、為線段長度；為從第e號(hào)三角形通過兩線段和流入內(nèi)的水量。根據(jù)（式3.1）的計(jì)算方法，對(duì)格點(diǎn)i周圍的所有三角形作類似計(jì)算，并求和得：（3.2） 另外，由于本文所研究的計(jì)算區(qū)域?yàn)槎ㄋ^邊界，承壓水，不考慮垂向補(bǔ)給量。 所以在單位時(shí)間內(nèi)，流入均衡區(qū)內(nèi)的總水量為：Q側(cè)。由水均衡原理得均衡方程：（i=1,2,…N） （3.3）式中：為的面積，為給水度。 這就是多邊形網(wǎng)格的水均衡方程。 由以上公式可以看出，在建立方程時(shí)，必須對(duì)格點(diǎn)i周圍的每個(gè)三角形計(jì)算、值，而在通常情況下，只知道格點(diǎn)的坐標(biāo)，因此，為了方便，下面將、以及均用格點(diǎn)坐標(biāo)標(biāo)出。（3.4）（3.5）并且記（3.6）（3.7）其中，=（3.8） 接下來計(jì)算多邊形的面積，可以將的計(jì)算分解成若干小三角形，分別計(jì)算每個(gè)小三角形的面積。計(jì)算三角形和三角形的面積。（3.9）同理可得（3.10） 類似此法，在格點(diǎn)i周圍的各個(gè)三角形中重復(fù)上述工作，即可求得的值。 由上述推導(dǎo)過程可見，要建立格點(diǎn)i的差分方程，只須對(duì)格點(diǎn)i周圍的三角形逐一做類似的計(jì)算即可。3.6.2三角形網(wǎng)格差分方程的解法 由（式3.2）和（式3.3）可知，在一個(gè)差分方程中，將涉及到格點(diǎn)i及其周圍幾個(gè)格點(diǎn)的水頭值。在本文的算例中，由于一類邊界格點(diǎn)是定水頭邊界，所以水頭值是已知的。其余內(nèi)格點(diǎn)是未知的，所以必須聯(lián)立方程組求解。 按前述方法，對(duì)每個(gè)內(nèi)格點(diǎn)建立差分方程，若內(nèi)格點(diǎn)有N1個(gè)，則共有N1個(gè)方程，而這些方程包含的未知數(shù)也正好是N1個(gè)，所以方程數(shù)和未知水頭數(shù)一樣多。此外，如果把格點(diǎn)i的差分方程中的項(xiàng)的系數(shù)排列在系數(shù)矩陣的主對(duì)角線上，其系數(shù)矩陣為對(duì)角占優(yōu)，故方程組存在唯一解，可用迭代法求解。 SOR（SuccessiveOver-Relaxation）法，即超松弛迭代法[22]，是目前解大型線型方程組的一種最常用的迭代加速方法[23]。它是在高斯-塞德爾迭代法[24]的基礎(chǔ)上進(jìn)行加速的，將前一步的結(jié)果與高斯-賽德爾迭代法的迭代值適當(dāng)?shù)募訖?quán)平均，期望獲得更好的近似值。其迭代公式[25][26]如下： i=1,2,……,n;k=0,1,2,… (3.11) SOR法中ω[27]的取值對(duì)迭代公式的收斂速度影響很大[28]，它的好壞直接影響到加速的快慢。為了保證迭代過程的收斂，必須要求0<ω<2，超松弛迭代法取1<ω<2[29]。算法3.4Step1： 將內(nèi)格點(diǎn)原始編號(hào)存入一個(gè)新的數(shù)組，獲得系數(shù)矩陣大小。Step2： 搜索三角形數(shù)組，找到包含當(dāng)前內(nèi)格點(diǎn)的三角形，計(jì)算當(dāng)前內(nèi)格點(diǎn)的多邊形區(qū)域面積，周圍格點(diǎn)的系數(shù)，存入數(shù)組中，直到所有內(nèi)格點(diǎn)都處理完。Step3： 將數(shù)組中的系數(shù)數(shù)據(jù)按照內(nèi)格點(diǎn)順序存入系數(shù)矩陣中，如果格點(diǎn)為邊界格點(diǎn)，將其計(jì)算結(jié)果放入右端已知項(xiàng)矩陣中。Step4： 使用SOR法解線性方程組，得到水頭值，存入數(shù)組中。3.6.3計(jì)算水頭模塊 三角形網(wǎng)格差分方程的解法在程序的“計(jì)算水頭”模塊中實(shí)現(xiàn)。該模塊前半部分為網(wǎng)格節(jié)點(diǎn)編號(hào)的優(yōu)化。后半部分為三角形網(wǎng)格有限差分法的計(jì)算。該部分主要的程序代碼為系數(shù)矩陣的形成。由于研究區(qū)域?yàn)槎ㄋ^邊界，水頭值已知，且不隨時(shí)間的變化而變化。所以這部分格點(diǎn)應(yīng)作為已知量，放到線性方程組右端已知項(xiàng)中。因此，用于計(jì)算的未知水頭值個(gè)數(shù)即為內(nèi)格點(diǎn)的個(gè)數(shù)。在已經(jīng)優(yōu)化的節(jié)點(diǎn)編號(hào)中需要除去邊界格點(diǎn)。僅保留內(nèi)格點(diǎn)。重新計(jì)算過的節(jié)點(diǎn)數(shù)組大小將只是內(nèi)格點(diǎn)的個(gè)數(shù)。以此數(shù)組的大小來定義系數(shù)矩陣，右端已知項(xiàng)矩陣，初始水頭值數(shù)組的大小。計(jì)算中包含的重要數(shù)據(jù)結(jié)構(gòu)包括：Matrix（points，points），即系數(shù)矩陣，BMatrix（points），即右端已知項(xiàng)矩陣，x（points），即初始水頭值向量，其中points即為內(nèi)格點(diǎn)的個(gè)數(shù)。以下描述程序的具體實(shí)現(xiàn)過程。首先形成新的節(jié)點(diǎn)數(shù)組，該數(shù)組不包含邊界格點(diǎn)，僅包含內(nèi)格點(diǎn)。這樣，也即確定了系數(shù)矩陣的大小，重新定義系數(shù)矩陣的大小，方程右端已知項(xiàng)的大小。搜索全局三角形數(shù)組，找到包含當(dāng)前節(jié)點(diǎn)的三角形，逐個(gè)計(jì)算當(dāng)前節(jié)點(diǎn)周圍節(jié)點(diǎn)的系數(shù)，存入當(dāng)前節(jié)點(diǎn)的系數(shù)數(shù)組中，計(jì)算當(dāng)前節(jié)點(diǎn)的多邊形區(qū)域面積，保存在變量中，待三角形全部搜索完畢，即獲得當(dāng)前節(jié)點(diǎn)多邊形區(qū)域的面積和周圍各個(gè)節(jié)點(diǎn)的系數(shù)。將數(shù)組中的各個(gè)節(jié)點(diǎn)的系數(shù)按照節(jié)點(diǎn)編號(hào)順序存入系數(shù)矩陣對(duì)應(yīng)位置，如果該節(jié)點(diǎn)為邊界節(jié)點(diǎn)，便獲得其水頭值，將水頭值與系數(shù)相乘加到右端已知項(xiàng)中。隨后將已知項(xiàng)存入已知項(xiàng)矩陣對(duì)應(yīng)位置。循環(huán)計(jì)算內(nèi)格點(diǎn)的系數(shù)和右端已知項(xiàng)，存入系數(shù)矩陣和右端已知項(xiàng)矩陣對(duì)應(yīng)位置，這樣便完成了系數(shù)矩陣和右端已知項(xiàng)矩陣的形成。取得各個(gè)格點(diǎn)的初始水頭值，存入初始向量數(shù)組中，將其和系數(shù)矩陣以及右端已知項(xiàng)矩陣作為SOR法的輸入數(shù)據(jù)。調(diào)用SOR法解方程組，迭代完成后即得到該時(shí)段末各個(gè)內(nèi)格點(diǎn)的水頭值。將其存入相應(yīng)數(shù)組中。以上，便是按照水均衡法原理[30]設(shè)計(jì)的礦井涌水量程序，程序獲得水頭值以后便可以交給相關(guān)部門預(yù)測涌水量，以及進(jìn)行礦井水害防治等工作。6 10-108工作面涌水量計(jì)算本文以霍州礦區(qū)白龍煤礦480水平，10-108工作面作為研究區(qū)域。將區(qū)域做了簡化，具體研究區(qū)域見圖6.1。并且假定工作面已經(jīng)形成，然后開始涌水，工作面不隨時(shí)間變化。圖6.1研究區(qū)域邊界及工作面位置外部多邊形區(qū)域?yàn)楸疚难芯繀^(qū)域的邊界，邊界為定水頭邊界，即水頭不隨時(shí)間變化，為一定值550。含水層為K2灰?guī)r，厚度為10m，給水度1/10000，滲透系數(shù)K為2m/天。內(nèi)部矩形區(qū)域?yàn)楹畬酉虏肯锏赖奈恢?，含水層中的水即?huì)從矩形區(qū)域向下涌入巷道，該區(qū)域內(nèi)部也為一定水頭邊界，為一定值480。滲透系數(shù)K為5m/天，下滲路徑長度為5m，矩形區(qū)域面積為100144m2。涌水量計(jì)算這也便是本文所要研究的主要內(nèi)容。首先對(duì)研究區(qū)域進(jìn)行初始化及剖分。結(jié)果如圖6.2、6.3所示。圖6.2研究區(qū)初始劃分圖6.3研究區(qū)Delaunay剖分 隨后建立差分方程組計(jì)算各個(gè)內(nèi)部格點(diǎn)的水頭值，然后應(yīng)用達(dá)西定律得到涌水量。達(dá)西定律如下： （6.1）式中：Q為涌水量，K為滲透系數(shù)，F(xiàn)為過水?dāng)嗝?，t為時(shí)間，h為總水頭損失，L為滲流路徑長度。 計(jì)算中得到的一些數(shù)據(jù)：多邊形區(qū)域邊界點(diǎn)個(gè)數(shù)：12工作面面積：100144m2剖分節(jié)點(diǎn)總個(gè)數(shù)：661三角形單元個(gè)數(shù)：1236內(nèi)部節(jié)點(diǎn)個(gè)數(shù)：649工作面內(nèi)節(jié)點(diǎn)個(gè)數(shù)：85 應(yīng)用程序得到的涌水量曲線如圖6.4，x軸為時(shí)間