王志東， 吳賀賀， 楊 爽, 陳明明

(江蘇科技大學 船舶與海洋工程學院, 江蘇 鎮(zhèn)江 212003)

動網(wǎng)格技術是數(shù)值分析大位移運動狀態(tài)下的物體流體力學性能的關鍵[1],常用的動網(wǎng)格技術包括網(wǎng)格移動與網(wǎng)格重新生成.網(wǎng)格移動有標準及改進的彈簧法和彈性體法,網(wǎng)格重新生成主要指網(wǎng)格的局部重構．

對于物體大位移運動的網(wǎng)格重構問題,重構生成網(wǎng)格的方法有Delaunay三角化方法和陣面推進法[2]．這兩種方法各有優(yōu)缺點,Delaunay三角化方法比陣面推進法具有較高的生成效率,同時生成的非結構網(wǎng)格接近正三角形的程度更高．但該方法比陣面推進法的布點能力弱,而布點能力的強弱是控制網(wǎng)格質(zhì)量的關鍵技術之一,所以當前的重構生成網(wǎng)格方法大都采用陣面推進法進行網(wǎng)格的重新生成[3]．文中提出一種改進方法，將Delaunay三角化方法應用到運動網(wǎng)格的重構中,充分利用Delaunay三角化方法的快速生成效率及接近正三角形的性質(zhì),同時提高Delaunay三角化方法的布點能力．

1 網(wǎng)格重構中Delaunay三角化改進方法

在進行網(wǎng)格重構時,一般在重構區(qū)域插入相同的單元數(shù),但對于物面變形較大的情況可能導致網(wǎng)格質(zhì)量較差[4].為保證網(wǎng)格質(zhì)量,文中以無量綱半徑作為網(wǎng)格生成的判斷依據(jù).無量綱半徑的引入可以充分提高網(wǎng)格質(zhì)量,并且可以提高網(wǎng)格生成程序的智能性．

1.1 三角形外接圓圓心長度標尺

(1)

式中：L(1),L(4),L(5)為邊界網(wǎng)格點的長度標尺;l1,l4,l5為Q點到3個頂點的距離．

1.2 三角形外接圓無量綱半徑

設任何一個Delaunay三角形為△k,外接圓半徑為r(k),外接圓圓心的長度標尺為l(k),則外接圓無量綱半徑R(k)為

(2)

圖1 初始化三角形R(r)?1的圖標

1.3 網(wǎng)格重構生成程序的主要步驟

PROGRAM MAIN

CALL creatgrid∥全場網(wǎng)格生成

CALL Reconstruction BoundAry∥確定重構邊界

DO t=1,n∥n為物面變形控制步長

CALL the deformAtion function of hill plane∥物面變形函數(shù)

CALL grid remeshing∥網(wǎng)格重構

CALL grid comBine∥網(wǎng)格合并

END DO

END

網(wǎng)格生成程序:

PROGRAM MAIN

CALL INPUT∥給定的初始條件

CALL INITIAL∥網(wǎng)格初始化

CALL NEWPOINT∥自動加點程序

CALL OUTPUT3∥網(wǎng)格信息輸出結果

END

文中在自動加點程序中將無量綱三角形半徑作為判斷是否繼續(xù)增加點的條件,保證生成的三角形最為接近正三角形．

網(wǎng)格生成中自動加點程序:

ENTRY NEWPOINT

DO WHILE(B(1).GT.S)∥S為設定的值

首先，農(nóng)業(yè)龍頭企業(yè)生產(chǎn)生活條件差，特別是種養(yǎng)企業(yè)生產(chǎn)環(huán)境差，沒有穩(wěn)定的職稱、收入與福利保障，難以吸引優(yōu)秀人才。

CALL CIRCLE∥求解三角形外接圓圓心

CALL SCAN∥判斷三角形外接圓圓心是否在計算區(qū)域內(nèi)

CALL INSERT∥插點程序

CALL SORTING∥插點后重新計算所有三角形外接圓無量綱半徑

END DO

END

在數(shù)組B(I)中存放的是所有三角形的無量綱半徑,且為由大到小的排列順序,則B(1)為最大的無量綱半徑．當B(1)大于設定值時自動加點,反之加點程序停止．

2 窗口邊界的確定方法

窗口邊界是指在動邊界周圍自動提取的一個封閉窗口,位于窗口邊界上和窗口邊界以外的網(wǎng)格點固定不動,而窗口內(nèi)部的網(wǎng)格可以變形或被重新生成．窗口邊界不僅是區(qū)分動點與不動點的分界線,也是采用陣面推進法重新生成變形區(qū)網(wǎng)格的初始陣面．因此,窗口邊界必須由現(xiàn)有網(wǎng)格單元的邊逐個連接構成,并且能夠在計算中自動、高效地提取并記錄．

窗口邊界的確定一般分為3種：① 由人為給定的到動邊界的法向距離確定;② 由包圍動邊界的基本幾何形狀確定;③ 由動邊界周圍的網(wǎng)格層數(shù)確定．其中第②、③種方法較為普遍．

為節(jié)省尋找窗口邊界的時間,在第②種方法的基礎上進行了改進.由于Delaunay三角化方法的布點能力比陣面推進法弱,為了提高其布點能力需要對網(wǎng)格進行加密處理,本文利用線元加密技術,同時將該加密邊界作為窗口邊界,這樣的改進可以避免復雜的尋點技術,而且利用常規(guī)尋點找出的邊界是鋸齒形狀的,這樣的復雜形狀在其邊界附近很難保證網(wǎng)格的密度．

現(xiàn)將兩種尋點算法歸納如下:

1) 常規(guī)算法

① 建立鏈表記錄每一個網(wǎng)格點的鄰點;

② 找到起始點L0,令START-POINT=LO;

③ 在START-POINT正方向查找最接近幾何輪廓線的鄰點L;

④ IF: L=L0 THEN

STOP

ELSE

令START-POINT=L,返回③

END IF

⑤ 將所得到的點加入到一維數(shù)組②中

2) 改進算法

改進算法中直接將已經(jīng)排列好的加密邊界點作為初始邊界,省去了進行尋點的時間．

3 用Delaunay三角化方法進行網(wǎng)格重構

網(wǎng)格重構的基本步驟:① Delaunay三角化方法在全場生成非結構網(wǎng)格;② 確定局部區(qū)域,確定方法主要是開“窗口”;③ 以“窗口”為邊界在內(nèi)部區(qū)域生成非結構網(wǎng)格,并求出窗口外部剩余網(wǎng)格;④ 將“窗口”內(nèi)網(wǎng)格與剩余網(wǎng)格合并．

在進行網(wǎng)格重構時,窗口的確定是關鍵步驟之一,文中在選擇“窗口”階段進行改進提高了網(wǎng)格質(zhì)量．此外,為了提高Delaunay三角化方法的布點能力,利用線元加密方法,在全場生成非結構網(wǎng)格時對網(wǎng)格進行加密,稱之為加密邊界,而該加密邊界作為重構“窗口”,不僅避免了復雜“窗口”的搜尋工作,而且增強了Delaunay三角化方法的布點能力,使Delaunay三角化方法在網(wǎng)格重構技術中得到應用．

4 網(wǎng)格質(zhì)量控制參數(shù)

4.1 網(wǎng)格長寬比和網(wǎng)格傾斜度

為進一步考察網(wǎng)格質(zhì)量,以網(wǎng)格長寬比和傾斜度為標準衡量網(wǎng)格質(zhì)量,其中網(wǎng)格長寬比定義為三角形網(wǎng)格外接圓半徑與內(nèi)切圓直徑之比,網(wǎng)格傾斜度定義為三角形網(wǎng)格內(nèi)角與等邊三角形內(nèi)角最大比值.

(3)

(4)

式中:Rout為外接圓半徑;Din為內(nèi)切圓半徑;Qc為等邊三角形內(nèi)角;即Qc=60°,Qmax,Qmin分別為三角形的最大內(nèi)角和最小內(nèi)角．由上式可知,網(wǎng)格長寬比越接近1,傾斜度越接近0,網(wǎng)格質(zhì)量越好．

4.2 形狀系數(shù)β

形狀系數(shù)β表征的是網(wǎng)格生成后偏離正三角形形狀的程度．正三角形的形狀系數(shù)為0．通過形狀系數(shù)可以判斷每個網(wǎng)格的質(zhì)量,β越接近0，網(wǎng)格質(zhì)量越高．形狀系數(shù)按下式計算:

(5)

(6)

(7)

5 算例分析

以柔性物面變形[6]為例,編寫了以局部重構方式生成的動網(wǎng)格程序．為了考查兩種網(wǎng)格重構方式的優(yōu)缺點,選取的區(qū)域是相同的,兩者窗口及加密邊界的選取相同,并選取相同時間刻度進行比較．

5.1 物面變形網(wǎng)格圖

圖2，3是在物面變形情況下的運動網(wǎng)格生成算例．圖2是經(jīng)過改進的網(wǎng)格局部重構方法,以初始的加密邊界為重構邊界.圖3是以常規(guī)方法尋找到的窗口作為重構邊界的網(wǎng)格局部生成方法．在生成網(wǎng)格時兩種方式中無量綱半徑設定值是相同的,但在插點過程中發(fā)現(xiàn),改進方法插點能力比傳統(tǒng)方法插點弱,特別是在窗口與加密邊界之間的網(wǎng)格密度兩者相差明顯,而且容易出現(xiàn)網(wǎng)格失效使程序無法運行．(T表示物面變形時刻,單位為s)

a) 外部剩余網(wǎng)格

b) T=0

c) T=15 s

a) 外部剩余網(wǎng)格

b) T=0

c) T=15 s

5.2 物面變形網(wǎng)格質(zhì)量分析圖

圖4，5給出了兩種方法下不同時刻物面發(fā)生變化后網(wǎng)格形狀系數(shù)的分布情況,β為形狀系數(shù).可以看出，在兩種方法中網(wǎng)格的形狀系數(shù)分布相差不大,但是常規(guī)方法比改進方法有少數(shù)點分布離散現(xiàn)象,個別網(wǎng)格嚴重偏離正三角形．

a) T=0

b) T=15 s

a) T=0

b) T=15 s

圖6,7中a)圖給出了兩種加密邊界[7]與窗口之間網(wǎng)格圖，b)圖給出了兩種方法在a)圖情況下的網(wǎng)格面積密度分布圖.從圖中得出,改進方法網(wǎng)格分布比較均勻,主要分布在0.000 05～0.000 7之間,而常規(guī)方法則出現(xiàn)了較明顯的偏離現(xiàn)象,最大值與最小值有成倍變化,說明改進方法在該區(qū)域網(wǎng)格密度變化較大．

a) 加密邊界與窗口之間網(wǎng)格

b) 加密與窗口之間網(wǎng)格面積函數(shù)分布

a) 加密邊界與窗口之間網(wǎng)格

b) 加密與窗口之間網(wǎng)格面積函數(shù)分布

圖8給出了兩種方法在不同時刻下網(wǎng)格長寬比和傾斜度的最大值、最小值以及平均值.可以看出在物面發(fā)生變化后網(wǎng)格質(zhì)量變化不大,兩種方法的最大和最小長寬比及傾斜度基本一致,改進方法能夠達到常規(guī)方法的精度．

圖8 兩種方法下網(wǎng)格在不同時間步長時的長寬比與傾斜度

6 結論

利用FORTRAN語言編制了兩種“窗口”邊界方式的局部重構法動網(wǎng)格生成程序,通過分析可得以下結論:

1) 利用加密邊界作為窗口邊界進行網(wǎng)格局部重構,不僅可以節(jié)省計算時間,而且增強了Delaunay三角化方法的布點能力;

2) 在相同物面變形條件下網(wǎng)格質(zhì)量得到了保證,有效增加了窗口與加密邊界之間的網(wǎng)格密度．

[1] Liu X Q, Li Q, Chai J Z.A new dynamic grid algrithm and its application [J].ACTAAeronauticaETAstronauticaSinica, 2008,29(4):817-822.

[2] 諸江. 非結構動網(wǎng)格生成方法研究[D].南京：南京理工大學,2006：2-3.

[3] 郭正. 包含運動邊界的多體非定常流場數(shù)值模擬方法研究 [D]. 湖南長沙：國防科學技術大學，2002：5-6.

[4] Liu Z S,Zhang Y F.Analysis and enhancement on lineal spring and torsional spring based dynamic mesh method[J].JournalofHarbinInstituteofTechnology,2005,37(8):1098-1102.

[5] 朱仁慶,李晏承,倪永燕,等.氣泡在水中上升運動的數(shù)值模擬[J].江蘇科技大學學報：自然科學版,2010,24(5):417-422.

Zhu Renqing,Li Yancheng,Ni Yongyan,et al.Numerical simulation of bubble rising in the water[J].JournalofJiangsuUniversityofScienceandTechnology:NaturalScienceEdition,2010,24(5):417-422.(in Chinese)

[6] 王志東,叢文超,李力軍.二維波狀擺動式魚類自主航行推進性能研究[J].江蘇科技大學學報:自然科學版,2010,24(3):217-222.

Wang Zhidong,Cong Wenchao,Li Lijun.Research on propulsion performance of autonomous navigation of 2D fish under the undulatory mode[J].JournalofJiangsuUniversityofScienceandTechnology:NaturalScienceEdition,2010,24(3):217-222.(in Chinese)

[7] Lohne R. Progress in grid generation via the advancing front technique[J].EngineeringwithComputers,1996, 12:186-210.