李小娟,孫建華
(江西省新干縣黃泥埠水庫管理局,江西 新干 331300)
隨著數(shù)值模擬技術在工程中的廣泛應用,人們對模擬計算精度的要求越來越高。網(wǎng)格的劃分是直接影響數(shù)值計算精度的一個重要因素。多年來,許多學者致力于網(wǎng)格的研究,并取得了豐碩的成果。計算網(wǎng)格通常劃分為結構網(wǎng)格和非結構網(wǎng)格。
結構網(wǎng)格包括代數(shù)方法中的多面法、無限插值法和通過求解等3種類型的微分方程生成的網(wǎng)格技術,結構網(wǎng)格具有編程快速、網(wǎng)格節(jié)點編號清晰等優(yōu)點,但當計算區(qū)域為天然河道時,邊界比較復雜,結構網(wǎng)格難以處理,生成網(wǎng)格的質(zhì)量難以把握。
非結構網(wǎng)格自20世紀80年代以來得到了迅猛的發(fā)展,Delaunay三角形剖分方法[1]和前沿推進法(advancing front method)在工程中得到了廣泛的應用。該方法對復雜邊界有很強的適應性,能生成質(zhì)量很高的三角形,但非結構網(wǎng)格點與點之間的編號十分雜亂,無規(guī)律可尋,難以在水沙二相流計算的離散格式中運用迎風性,不利于水流的模擬計算[2-3]。
本文充分發(fā)揮了結構網(wǎng)格和非結構網(wǎng)格的優(yōu)點,提出一種新的網(wǎng)格生成方法。該方法采用使網(wǎng)格沿水流的方向由上游向下游逐步生成,在網(wǎng)格的自動生成過程中以前沿推進法的思想為主導,引入Delaunay三角形的概念,用以控制三角形的質(zhì)量,最后生成優(yōu)質(zhì)的三角形網(wǎng)格。
本河道水流模擬研究中,取一進口斷面和一出口斷面,加上河道的天然邊界,形成一封閉的區(qū)域,圖1所示。
圖1 河道邊界離散示意圖
圖1中,AB為進口斷面,CD為出口斷面,BC和AD為河道的天然邊界,要在ABCD區(qū)域中劃分網(wǎng)格,首先要對ABCD的邊界進行離散。分別給出A、B、C、D的坐標(XA,YA),(XB,YB),(XC,YC),(XD,YD),引入網(wǎng)格尺度[4]。網(wǎng)格尺度為某一段要求劃分網(wǎng)格的尺寸。設A、B、C、D 的網(wǎng)格尺度分別為DA、DB、DC、DD,設定了網(wǎng)格尺度標志著邊界離散的尺寸已經(jīng)確定,邊界的點數(shù)和坐標可計算出來。
區(qū)域內(nèi)節(jié)點和單元的生成如下:
(1)將AB邊作為進口邊,AB上離散點為前沿點,由相鄰前沿點所組成的邊為前沿邊,從前沿邊向下游擴展。
(2)選定 AB邊上的第一條離散邊 A2,其中“2”為頂點名稱,其邊長為LA2,在有向線段中點向右邊作垂直于的線段,線段長,H 點為擴展的一新點,其示意圖如圖2。
圖2 前沿邊擴展示意圖
(3)判斷點H是否為符合要求的擴展點。以H為圓心,以 αLAM為半徑劃圓(α 一般取 0.5~0.8),判斷是否有離散點或前沿點被包括在圓內(nèi),若圓內(nèi)沒有前沿點,則(2)中所得的 H 點為所擴展的點,H 為新的前沿點,AH、H2為新的前沿邊。如果圓周內(nèi)有離散點,如圖3。求出除H外的圓內(nèi)各點對的望角,即∠AI2、∠AJ2、∠AK2、∠AL2、∠AM2,比較各望角的大小,顯然∠AI2最大,于是選I為所擴展的新點,如果在 BC上,新增的前沿邊為;如果在DA上,新增的前沿邊為;如果點 I在上,新增的前沿邊為和。至此邊的擴展完成,重復步驟(1)~(3)對邊進行擴展,直到最后一個前沿邊擴展完畢。
(4)刪除編號重合的三角形。在上述擴展過程中同一個前沿邊有可能擴展了多次,也就是說對于同一個三角形可能重復計算了。比較各單元的節(jié)點編號,如果幾個單元的節(jié)點編號完全相同則僅留下單元號最小的單元,其余的單元均刪除。如45單元的節(jié)點編號分別為22、23、30;50 單元的節(jié)點編號分別為30、22、23。 兩單元的編號完全相同,由于45<50,故刪掉50單元,只保留45單元,其余單元數(shù)大于50的單元其單元數(shù)均減1。如此類推,最終區(qū)域ABCD內(nèi)沒有重合的單元。
圖3 確定最終擴展點示意圖
三角形網(wǎng)格優(yōu)化的主要目的是使網(wǎng)格接近正三角形,保證網(wǎng)格的質(zhì)量。網(wǎng)格的優(yōu)化方法主要有“結構優(yōu)化”和“位置優(yōu)化”。前者調(diào)整網(wǎng)的拓撲結構,后者調(diào)整內(nèi)部節(jié)點的位置。
在文獻[5]中作者引入了節(jié)點的“度”的概念,定義度為共享該節(jié)點的單元數(shù)目。并得出三角形節(jié)點的度δ一般小于8,最佳取值范圍為5≤δ≤8。因此,僅需要對δ等與3或4的三角形節(jié)點進行處理。
判斷三角形節(jié)點的度,邊界點除外,當δ=3時,去掉該節(jié)點,其余三邊拼成一新的三角形,如圖4。去掉點D,ABC組成新的三角形。
圖4 δ=3 時轉換圖
當δ=4時計算各單元中該角的角度,當角度大于100時需要添加邊。如圖5,假設∠AEB>100;∠BEC>100,在∠AEB和∠BEC間要添加一邊,AB,BC分別為∠AEB和∠BEC對應的邊,取出共享AB、BC邊的三角形ΔAGB和ΔBFC,分別消去邊AB和BC,連接EG和EF,使三角形得以優(yōu)化。
圖5 δ=4 時轉換圖
網(wǎng)格的位置優(yōu)化通常用Laplacian均勻計算公式[6]:
XN、YN—分別為與I點連接的各點的橫坐標和縱坐標。
上述公式是簡單的線性平均,使I點處于所在面的幾何中心。該計算公式在運用過程中變化很大,很有可能移到某個相鄰的單元內(nèi),并且完全忽視了原來(XI,YI)對它的影響,為此對該種方法進行了改進,引入松弛迭代技術:
式中:(XI,YI)—原來 I點的坐標;
NI—與I點相連的點的個數(shù),如圖5中,與D點相連的點為A、B、C共有3個;
(XN,YN)—相連點的坐標;
ω—松弛因子,可以取1。
對于網(wǎng)格位置的優(yōu)化可以進行多次光滑,每次光滑均取最新坐標,一般經(jīng)過5~10次可以達到較好的效果。
水流連續(xù)方程:
式中:U、V—分別為垂線平均流速在x、y方向上的分量;
ZS、Zb—分別為水位和河底高程;
H—垂線水深;
ρ—水密度;
νt—紊流粘滯性系數(shù);
τx、τy—分別為底部切應力在 x、y方向上的分量;
f—柯氏力系數(shù),f=2WsinΦ;(W 為地球自轉角速度,Φ為地理緯度)。
其中,
式中:C—謝才系數(shù)。
對于不動邊界條件,可分為以下4類處理:
(1)上游進口邊界(開邊界)Γ1,給出
(2)下游出口邊界(開邊界)Γ2,給出
(3)不滑動岸壁邊界(閉邊界)Γ3,給出
(4)滑動岸壁邊界(閉邊界)Γ4,規(guī)定
如圖6為概化的平面河道,河道長2200 m,寬600 m,河道中存在一座小島,河底比降設定為0.0001。采用本文提出的三角形網(wǎng)格自動生成方法,該區(qū)域共生成684個節(jié)點,1236個網(wǎng)格單元,網(wǎng)格長約50 m。當進口斷面給定流量3000 m3/s,下游給定恒定水位3 m時進行計算。圖7為根據(jù)生成的網(wǎng)格計算的流場圖,可以看出水流從進口斷面向下游流動,遇小島后水流分流,在小島后形成回流,且流速較小,該流速變化規(guī)律與實際情況相符。
圖6 河道網(wǎng)格劃分
圖7 河道平面流場圖
本文將波前推進法和Delaunay三角形剖分方法結合起來,選取進口斷面為網(wǎng)格擴展的前沿,提出了適用于平面二維河道水流模擬的網(wǎng)格自動剖分方法,在平面河道中生成了優(yōu)質(zhì)的三角形網(wǎng)格。該方法能快速實現(xiàn)河道網(wǎng)格的自動剖分,生成后的網(wǎng)格具有優(yōu)質(zhì)、對邊界的適應性強等特點,能提供數(shù)值模擬計算需要的前處理數(shù)據(jù)?;谏傻木W(wǎng)格應用于水流模擬計算時,計算結果能充分反映水流變化的實際情況。
[1]楊曉東,劉春太,申長雨,等.內(nèi)部帶特征約束的任意平面域的三角形網(wǎng)格生成方法[J].計算物理,2000(3).
[2]張細兵.河道有限元網(wǎng)格自動剖分方法的研究[J].長江科學院院報,2000(3).
[3]吳 騰,鐘德鈺,張紅武.水庫自適應控制運用模型及其在亭口水庫的應用[J].水力發(fā)電學報,2010,29(3):97-102.
[4]張 征,李孟國.三角形網(wǎng)格自動生成技術的應用[J].水道港口,2001(3).
[5]張均鋒,劉桂齋,陳 剛.有限元平面元三角形網(wǎng)格的優(yōu)化[J].山東礦業(yè)學院學學報,1997(3).
[6]羅特軍,羅季軍,汪 榴,等.有限元網(wǎng)格優(yōu)化方法[J],四川聯(lián)合大學學報,1999(3).