三角網(wǎng)格剖分方法及其在河道水流模擬中的應(yīng)用

李小娟，孫建華

（江西省新干縣黃泥埠水庫管理局，江西 新干 331300）

0 引 言

隨著數(shù)值模擬技術(shù)在工程中的廣泛應(yīng)用，人們對模擬計(jì)算精度的要求越來越高。網(wǎng)格的劃分是直接影響數(shù)值計(jì)算精度的一個重要因素。多年來，許多學(xué)者致力于網(wǎng)格的研究，并取得了豐碩的成果。計(jì)算網(wǎng)格通常劃分為結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格。

結(jié)構(gòu)網(wǎng)格包括代數(shù)方法中的多面法、無限插值法和通過求解等3種類型的微分方程生成的網(wǎng)格技術(shù)，結(jié)構(gòu)網(wǎng)格具有編程快速、網(wǎng)格節(jié)點(diǎn)編號清晰等優(yōu)點(diǎn)，但當(dāng)計(jì)算區(qū)域?yàn)樘烊缓拥罆r，邊界比較復(fù)雜，結(jié)構(gòu)網(wǎng)格難以處理，生成網(wǎng)格的質(zhì)量難以把握。

非結(jié)構(gòu)網(wǎng)格自20世紀(jì)80年代以來得到了迅猛的發(fā)展，Delaunay三角形剖分方法[1]和前沿推進(jìn)法（advancing front method）在工程中得到了廣泛的應(yīng)用。該方法對復(fù)雜邊界有很強(qiáng)的適應(yīng)性，能生成質(zhì)量很高的三角形，但非結(jié)構(gòu)網(wǎng)格點(diǎn)與點(diǎn)之間的編號十分雜亂，無規(guī)律可尋，難以在水沙二相流計(jì)算的離散格式中運(yùn)用迎風(fēng)性，不利于水流的模擬計(jì)算[2-3]。

本文充分發(fā)揮了結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格的優(yōu)點(diǎn)，提出一種新的網(wǎng)格生成方法。該方法采用使網(wǎng)格沿水流的方向由上游向下游逐步生成，在網(wǎng)格的自動生成過程中以前沿推進(jìn)法的思想為主導(dǎo)，引入Delaunay三角形的概念，用以控制三角形的質(zhì)量，最后生成優(yōu)質(zhì)的三角形網(wǎng)格。

1 河道邊界的處理

本河道水流模擬研究中，取一進(jìn)口斷面和一出口斷面，加上河道的天然邊界，形成一封閉的區(qū)域，圖1所示。

圖1 河道邊界離散示意圖

圖1中，AB為進(jìn)口斷面，CD為出口斷面，BC和AD為河道的天然邊界，要在ABCD區(qū)域中劃分網(wǎng)格，首先要對ABCD的邊界進(jìn)行離散。分別給出A、B、C、D的坐標(biāo)（XA，YA），（XB，YB），（XC，YC），（XD，YD），引入網(wǎng)格尺度[4]。網(wǎng)格尺度為某一段要求劃分網(wǎng)格的尺寸。設(shè)A、B、C、D 的網(wǎng)格尺度分別為DA、DB、DC、DD，設(shè)定了網(wǎng)格尺度標(biāo)志著邊界離散的尺寸已經(jīng)確定，邊界的點(diǎn)數(shù)和坐標(biāo)可計(jì)算出來。

2 區(qū)域內(nèi)節(jié)點(diǎn)和單元的生成

區(qū)域內(nèi)節(jié)點(diǎn)和單元的生成如下：

（1）將AB邊作為進(jìn)口邊，AB上離散點(diǎn)為前沿點(diǎn)，由相鄰前沿點(diǎn)所組成的邊為前沿邊，從前沿邊向下游擴(kuò)展。

（2）選定 AB邊上的第一條離散邊 A2，其中“2”為頂點(diǎn)名稱，其邊長為LA2，在有向線段中點(diǎn)向右邊作垂直于的線段，線段長，H 點(diǎn)為擴(kuò)展的一新點(diǎn)，其示意圖如圖2。

圖2 前沿邊擴(kuò)展示意圖

（3）判斷點(diǎn)H是否為符合要求的擴(kuò)展點(diǎn)。以H為圓心，以 αLAM為半徑劃圓（α 一般取 0.5～0.8），判斷是否有離散點(diǎn)或前沿點(diǎn)被包括在圓內(nèi)，若圓內(nèi)沒有前沿點(diǎn)，則（2）中所得的 H 點(diǎn)為所擴(kuò)展的點(diǎn)，H 為新的前沿點(diǎn)，AH、H2為新的前沿邊。如果圓周內(nèi)有離散點(diǎn)，如圖3。求出除H外的圓內(nèi)各點(diǎn)對的望角，即∠AI2、∠AJ2、∠AK2、∠AL2、∠AM2，比較各望角的大小，顯然∠AI2最大，于是選I為所擴(kuò)展的新點(diǎn)，如果在 BC上，新增的前沿邊為；如果在DA上，新增的前沿邊為；如果點(diǎn) I在上，新增的前沿邊為和。至此邊的擴(kuò)展完成，重復(fù)步驟（1）～（3）對邊進(jìn)行擴(kuò)展，直到最后一個前沿邊擴(kuò)展完畢。

（4）刪除編號重合的三角形。在上述擴(kuò)展過程中同一個前沿邊有可能擴(kuò)展了多次，也就是說對于同一個三角形可能重復(fù)計(jì)算了。比較各單元的節(jié)點(diǎn)編號，如果幾個單元的節(jié)點(diǎn)編號完全相同則僅留下單元號最小的單元，其余的單元均刪除。如45單元的節(jié)點(diǎn)編號分別為22、23、30；50 單元的節(jié)點(diǎn)編號分別為30、22、23。 兩單元的編號完全相同，由于45＜50，故刪掉50單元，只保留45單元，其余單元數(shù)大于50的單元其單元數(shù)均減1。如此類推，最終區(qū)域ABCD內(nèi)沒有重合的單元。

圖3 確定最終擴(kuò)展點(diǎn)示意圖

3 網(wǎng)格的光滑處理

三角形網(wǎng)格優(yōu)化的主要目的是使網(wǎng)格接近正三角形，保證網(wǎng)格的質(zhì)量。網(wǎng)格的優(yōu)化方法主要有“結(jié)構(gòu)優(yōu)化”和“位置優(yōu)化”。前者調(diào)整網(wǎng)的拓?fù)浣Y(jié)構(gòu)，后者調(diào)整內(nèi)部節(jié)點(diǎn)的位置。

3.1 網(wǎng)格的結(jié)構(gòu)優(yōu)化

在文獻(xiàn)[5]中作者引入了節(jié)點(diǎn)的“度”的概念，定義度為共享該節(jié)點(diǎn)的單元數(shù)目。并得出三角形節(jié)點(diǎn)的度δ一般小于8，最佳取值范圍為5≤δ≤8。因此，僅需要對δ等與3或4的三角形節(jié)點(diǎn)進(jìn)行處理。

判斷三角形節(jié)點(diǎn)的度，邊界點(diǎn)除外，當(dāng)δ＝3時，去掉該節(jié)點(diǎn)，其余三邊拼成一新的三角形，如圖4。去掉點(diǎn)D，ABC組成新的三角形。

圖4 δ＝3 時轉(zhuǎn)換圖

當(dāng)δ＝4時計(jì)算各單元中該角的角度，當(dāng)角度大于100時需要添加邊。如圖5，假設(shè)∠AEB＞100；∠BEC＞100，在∠AEB和∠BEC間要添加一邊，AB，BC分別為∠AEB和∠BEC對應(yīng)的邊，取出共享AB、BC邊的三角形ΔAGB和ΔBFC，分別消去邊AB和BC，連接EG和EF，使三角形得以優(yōu)化。

圖5 δ＝4 時轉(zhuǎn)換圖

3.2 網(wǎng)格的位置優(yōu)化

網(wǎng)格的位置優(yōu)化通常用Laplacian均勻計(jì)算公式[6]：

XN、YN—分別為與I點(diǎn)連接的各點(diǎn)的橫坐標(biāo)和縱坐標(biāo)。

上述公式是簡單的線性平均，使I點(diǎn)處于所在面的幾何中心。該計(jì)算公式在運(yùn)用過程中變化很大，很有可能移到某個相鄰的單元內(nèi)，并且完全忽視了原來（XI，YI）對它的影響，為此對該種方法進(jìn)行了改進(jìn)，引入松弛迭代技術(shù)：

式中：（XI，YI）—原來 I點(diǎn)的坐標(biāo)；

NI—與I點(diǎn)相連的點(diǎn)的個數(shù)，如圖5中，與D點(diǎn)相連的點(diǎn)為A、B、C共有3個；

（XN，YN）—相連點(diǎn)的坐標(biāo)；

ω—松弛因子，可以取1。

對于網(wǎng)格位置的優(yōu)化可以進(jìn)行多次光滑，每次光滑均取最新坐標(biāo)，一般經(jīng)過5～10次可以達(dá)到較好的效果。

4 三角網(wǎng)格剖分方法的應(yīng)用

4.1 基本方程

水流連續(xù)方程：

式中：U、V—分別為垂線平均流速在x、y方向上的分量；

ZS、Zb—分別為水位和河底高程；

H—垂線水深；

ρ—水密度；

νt—紊流粘滯性系數(shù)；

τx、τy—分別為底部切應(yīng)力在 x、y方向上的分量；

f—柯氏力系數(shù)，f＝2WsinΦ；（W 為地球自轉(zhuǎn)角速度，Φ為地理緯度）。

其中，

式中：C—謝才系數(shù)。

4.2 邊界條件處理

對于不動邊界條件，可分為以下4類處理：

（1）上游進(jìn)口邊界（開邊界）Γ1，給出

（2）下游出口邊界（開邊界）Γ2，給出

（3）不滑動岸壁邊界（閉邊界）Γ3，給出

（4）滑動岸壁邊界（閉邊界）Γ4，規(guī)定

4.3 計(jì)算實(shí)例

如圖6為概化的平面河道，河道長2200 m，寬600 m，河道中存在一座小島，河底比降設(shè)定為0.0001。采用本文提出的三角形網(wǎng)格自動生成方法，該區(qū)域共生成684個節(jié)點(diǎn)，1236個網(wǎng)格單元，網(wǎng)格長約50 m。當(dāng)進(jìn)口斷面給定流量3000 m3/s，下游給定恒定水位3 m時進(jìn)行計(jì)算。圖7為根據(jù)生成的網(wǎng)格計(jì)算的流場圖，可以看出水流從進(jìn)口斷面向下游流動，遇小島后水流分流，在小島后形成回流，且流速較小，該流速變化規(guī)律與實(shí)際情況相符。

圖6 河道網(wǎng)格劃分

圖7 河道平面流場圖

5 結(jié)語

本文將波前推進(jìn)法和Delaunay三角形剖分方法結(jié)合起來，選取進(jìn)口斷面為網(wǎng)格擴(kuò)展的前沿，提出了適用于平面二維河道水流模擬的網(wǎng)格自動剖分方法，在平面河道中生成了優(yōu)質(zhì)的三角形網(wǎng)格。該方法能快速實(shí)現(xiàn)河道網(wǎng)格的自動剖分，生成后的網(wǎng)格具有優(yōu)質(zhì)、對邊界的適應(yīng)性強(qiáng)等特點(diǎn)，能提供數(shù)值模擬計(jì)算需要的前處理數(shù)據(jù)?；谏傻木W(wǎng)格應(yīng)用于水流模擬計(jì)算時，計(jì)算結(jié)果能充分反映水流變化的實(shí)際情況。

［1］楊曉東，劉春太，申長雨，等.內(nèi)部帶特征約束的任意平面域的三角形網(wǎng)格生成方法［J］.計(jì)算物理，2000（3）.

［2］張細(xì)兵.河道有限元網(wǎng)格自動剖分方法的研究［J］.長江科學(xué)院院報，2000（3）.

［3］吳 騰，鐘德鈺，張紅武.水庫自適應(yīng)控制運(yùn)用模型及其在亭口水庫的應(yīng)用［J］.水力發(fā)電學(xué)報，2010，29（3）：97－102.

［4］張 征，李孟國.三角形網(wǎng)格自動生成技術(shù)的應(yīng)用［J］.水道港口，2001（3）.

［5］張均鋒，劉桂齋，陳 剛.有限元平面元三角形網(wǎng)格的優(yōu)化［J］.山東礦業(yè)學(xué)院學(xué)學(xué)報，1997（3）.

［6］羅特軍，羅季軍，汪 榴，等.有限元網(wǎng)格優(yōu)化方法［J］，四川聯(lián)合大學(xué)學(xué)報，1999（3）.