動(dòng)網(wǎng)格技術(shù)中局部重構(gòu)方法的改進(jìn)研究

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

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

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

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

1 網(wǎng)格重構(gòu)中Delaunay三角化改進(jìn)方法

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

1.1 三角形外接圓圓心長(zhǎng)度標(biāo)尺

(1)

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

1.2 三角形外接圓無(wú)量綱半徑

設(shè)任何一個(gè)Delaunay三角形為△k,外接圓半徑為r(k),外接圓圓心的長(zhǎng)度標(biāo)尺為l(k),則外接圓無(wú)量綱半徑R(k)為

(2)

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

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

PROGRAM MAIN

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

CALL Reconstruction BoundAry∥確定重構(gòu)邊界

DO t=1,n∥n為物面變形控制步長(zhǎng)

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

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

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

END DO

END

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

PROGRAM MAIN

CALL INPUT∥給定的初始條件

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

CALL NEWPOINT∥自動(dòng)加點(diǎn)程序

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

END

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

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

ENTRY NEWPOINT

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

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

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

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

CALL INSERT∥插點(diǎn)程序

CALL SORTING∥插點(diǎn)后重新計(jì)算所有三角形外接圓無(wú)量綱半徑

END DO

END

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

2 窗口邊界的確定方法

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

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

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

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

1) 常規(guī)算法

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

② 找到起始點(diǎn)L0,令START-POINT=LO;

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

④ IF: L=L0 THEN

STOP

ELSE

令START-POINT=L,返回③

END IF

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

2) 改進(jìn)算法

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

3 用Delaunay三角化方法進(jìn)行網(wǎng)格重構(gòu)

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

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

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

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

為進(jìn)一步考察網(wǎng)格質(zhì)量,以網(wǎng)格長(zhǎng)寬比和傾斜度為標(biāo)準(zhǔn)衡量網(wǎng)格質(zhì)量,其中網(wǎng)格長(zhǎ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)格長(zhǎng)寬比越接近1,傾斜度越接近0,網(wǎng)格質(zhì)量越好．

4.2 形狀系數(shù)β

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

(5)

(6)

(7)

5 算例分析

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

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

圖2，3是在物面變形情況下的運(yùn)動(dòng)網(wǎng)格生成算例．圖2是經(jīng)過(guò)改進(jìn)的網(wǎng)格局部重構(gòu)方法,以初始的加密邊界為重構(gòu)邊界.圖3是以常規(guī)方法尋找到的窗口作為重構(gòu)邊界的網(wǎng)格局部生成方法．在生成網(wǎng)格時(shí)兩種方式中無(wú)量綱半徑設(shè)定值是相同的,但在插點(diǎn)過(guò)程中發(fā)現(xiàn),改進(jìn)方法插點(diǎn)能力比傳統(tǒng)方法插點(diǎn)弱,特別是在窗口與加密邊界之間的網(wǎng)格密度兩者相差明顯,而且容易出現(xiàn)網(wǎng)格失效使程序無(wú)法運(yùn)行．(T表示物面變形時(shí)刻,單位為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給出了兩種方法下不同時(shí)刻物面發(fā)生變化后網(wǎng)格形狀系數(shù)的分布情況,β為形狀系數(shù).可以看出，在兩種方法中網(wǎng)格的形狀系數(shù)分布相差不大,但是常規(guī)方法比改進(jìn)方法有少數(shù)點(diǎn)分布離散現(xiàn)象,個(gè)別網(wǎng)格嚴(yán)重偏離正三角形．

a) T=0

b) T=15 s

a) T=0

b) T=15 s

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

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

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

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

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

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

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

6 結(jié)論

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

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

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] 諸江. 非結(jié)構(gòu)動(dòng)網(wǎng)格生成方法研究[D].南京：南京理工大學(xué),2006：2-3.

[3] 郭正. 包含運(yùn)動(dòng)邊界的多體非定常流場(chǎng)數(shù)值模擬方法研究 [D]. 湖南長(zhǎng)沙：國(guó)防科學(xué)技術(shù)大學(xué)，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] 朱仁慶,李晏承,倪永燕,等.氣泡在水中上升運(yùn)動(dòng)的數(shù)值模擬[J].江蘇科技大學(xué)學(xué)報(bào)：自然科學(xué)版,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] 王志東,叢文超,李力軍.二維波狀擺動(dòng)式魚(yú)類自主航行推進(jìn)性能研究[J].江蘇科技大學(xué)學(xué)報(bào):自然科學(xué)版,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.