基于直角切割點(diǎn)調(diào)整的二維網(wǎng)格生成方法

于上

摘要：目前網(wǎng)格生成占用了流場計(jì)算的大量時(shí)間，為了進(jìn)一步縮短網(wǎng)格生成的計(jì)算時(shí)間，該文在直角網(wǎng)格上直接進(jìn)行快速的網(wǎng)格調(diào)整，針對多種邊界情況進(jìn)行了分類，保持了對復(fù)雜邊界良好的適應(yīng)性，同時(shí)在網(wǎng)格質(zhì)量可以接受的范圍內(nèi)，可以在與網(wǎng)格單元數(shù)量成線性增長的時(shí)間復(fù)雜度下生成網(wǎng)格。最后以NACA0012翼型為例給出計(jì)算結(jié)果，驗(yàn)證了算法的有效性。

關(guān)鍵詞：網(wǎng)格生成；直角網(wǎng)格；網(wǎng)格邊界調(diào)整

中圖分類號：TP18 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號：1009-3044（2015）22-0164-03

Abstract： Currently， mesh generation for CFD takes up considerable time of the whole calculation. To further reduce the time spent on mesh generation， a new method has been proposed to make adjustment quickly on Cartesian mesh for complex shapes. The mesh is adjusted according to different situations of mesh boundary， so it performs well on complex boundary. The quality of mesh will be in an acceptable range and the time complexity is linear to the number of grid cells. An example of generation on NACA0012 airfoil is given to confirm the effectiveness of the algorithm.

Key words： mesh generation； Cartesian mesh； mesh boundary adjustment

網(wǎng)格生成技術(shù)，是計(jì)算流體力學(xué)（CFD）的關(guān)鍵技術(shù)之一。網(wǎng)格生成的過程，把連續(xù)的計(jì)算區(qū)域劃分成許多子區(qū)域，以便進(jìn)行后續(xù)的計(jì)算。在目前的CFD計(jì)算周期中，網(wǎng)格生成占用了很大一部分時(shí)間，并且后續(xù)的計(jì)算精度和計(jì)算效率也很大程度上依賴于網(wǎng)格生成的質(zhì)量和方法。

生成的網(wǎng)格從拓?fù)浣Y(jié)構(gòu)上分，主要分為結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格。結(jié)構(gòu)化網(wǎng)格數(shù)據(jù)結(jié)構(gòu)簡單，生成速度快。結(jié)構(gòu)化網(wǎng)格在網(wǎng)格邊界處通常使用適體坐標(biāo)法[1]，但對于復(fù)雜外形的適應(yīng)性比較差。非結(jié)構(gòu)化網(wǎng)格對于復(fù)雜區(qū)域的適應(yīng)性較強(qiáng)，但網(wǎng)格質(zhì)量不好控制，計(jì)算效率相對結(jié)構(gòu)化網(wǎng)格要低?，F(xiàn)在非結(jié)構(gòu)網(wǎng)格通常采用帶約束的Delaunay三角剖分的方法[2]或四叉樹方法[3]。針對不同網(wǎng)格的不同優(yōu)點(diǎn)及缺點(diǎn)，混合網(wǎng)格應(yīng)用的也越來越廣泛[4]。

為了進(jìn)一步提升網(wǎng)格生成的計(jì)算效率，本文所采用了一種在直角網(wǎng)格上調(diào)整頂點(diǎn)來適應(yīng)邊界的生成方法。這樣做的優(yōu)點(diǎn)是：與背景的結(jié)構(gòu)化網(wǎng)格銜接時(shí)不需要做特殊處理；在邊界附近對被切割的網(wǎng)格分情況討論去調(diào)整網(wǎng)格，利用直角網(wǎng)格生成速度快的優(yōu)勢，避免了復(fù)雜度較高的Delaunay三角剖分。

本文的方法主要分為3個(gè)步驟，生成邊界網(wǎng)格與計(jì)算交點(diǎn)，調(diào)整交點(diǎn)位置，生成最終網(wǎng)格，最后以NACA0012翼型為例，給出網(wǎng)格的生成結(jié)果。

1 光滑邊界的網(wǎng)格生成

Wenjun Ying和Wei-Cheng Wang的方法[5]可以直接應(yīng)用到光滑邊界的網(wǎng)格生成中。在該方法中，所有邊界與網(wǎng)格的交點(diǎn)都被稱為控制點(diǎn)。與網(wǎng)格點(diǎn)距離近的稱為主要控制點(diǎn)，與網(wǎng)格點(diǎn)距離遠(yuǎn)的稱為次要控制點(diǎn)。

將結(jié)構(gòu)網(wǎng)格上頂點(diǎn) 坐標(biāo)移動(dòng)到對應(yīng)的主要控制點(diǎn)，就完成了網(wǎng)格的調(diào)整，如圖1（b）的箭頭方向所示。

然而這個(gè)方法并沒有考慮到復(fù)雜外形邊界上不光滑的尖點(diǎn)。所以本文在此基礎(chǔ)上，針對尖點(diǎn)的邊界附近的所有情況給出調(diào)整方案，把這個(gè)調(diào)整邊界網(wǎng)格的方法推廣到復(fù)雜的外形上。

2 尖點(diǎn)附近的網(wǎng)格生成

這里的“尖點(diǎn)”是指曲線在該點(diǎn)處，兩側(cè)的導(dǎo)數(shù)不相等。

2.1 尖點(diǎn)附近需要處理的問題

光滑邊界的外形所有的點(diǎn)附近都是平滑的，而尖點(diǎn)兩側(cè)曲線的切向不同，如果僅僅用交點(diǎn)去近似，會(huì)非常不準(zhǔn)確。

假設(shè)網(wǎng)格足夠密，使得外形在3*3的網(wǎng)格范圍內(nèi)最多只有1個(gè)尖點(diǎn)。又因?yàn)檫@個(gè)頂點(diǎn)兩端的切向量不同，所以對于方形結(jié)構(gòu)網(wǎng)格的每個(gè)頂點(diǎn)，最多有2條近似的直線與頂點(diǎn)周圍的4條邊產(chǎn)生交點(diǎn)。

針對2.1中的頂點(diǎn)以及圖2的6種情況，將會(huì)在2.3分別給出調(diào)整的方案。在2.2中先給出了計(jì)算前的預(yù)處理步驟。

2.2 網(wǎng)格調(diào)整前的預(yù)處理

初始的網(wǎng)格是與背景網(wǎng)格相同的正方形網(wǎng)格，根據(jù)四個(gè)頂點(diǎn)在網(wǎng)格的內(nèi)外的情況，劃分出流場網(wǎng)格、物面網(wǎng)格和固體網(wǎng)格。計(jì)算區(qū)域在現(xiàn)有的物面網(wǎng)格的基礎(chǔ)上，向外側(cè)擴(kuò)展2層，擴(kuò)展后的物面網(wǎng)格的外側(cè)頂點(diǎn)會(huì)自然地與背景網(wǎng)格銜接起來。之后計(jì)算外形與網(wǎng)格線的交點(diǎn)。

2.3 直角網(wǎng)格的調(diào)整方法

2.3.1 尖點(diǎn)的位置調(diào)整

對于外形上的尖點(diǎn)，穿過尖點(diǎn)的曲線兩側(cè)切向不同，不可能用一條直線近似。所以調(diào)整的方式是直接把距離最近的網(wǎng)格點(diǎn)移動(dòng)到這個(gè)尖點(diǎn)處，如圖3所示。

對網(wǎng)格調(diào)整后，最終的生成結(jié)果還需要去除內(nèi)部的網(wǎng)格以及生成三角形網(wǎng)格兩步工作。

2.4 去除內(nèi)部的網(wǎng)格

如果外形沒有凹角，4個(gè)點(diǎn)都不在外部即可判為內(nèi)部網(wǎng)格。允許外形有凹角時(shí)，內(nèi)部網(wǎng)格的判定分為如下幾種：

第1種是4個(gè)頂點(diǎn)都不在外部，且有2個(gè)或2個(gè)以上是內(nèi)部點(diǎn)，一定是內(nèi)部網(wǎng)格。

第2種是1個(gè)頂點(diǎn)在內(nèi)部，3個(gè)頂點(diǎn)在邊界上的情況，此時(shí)如果網(wǎng)格的2條邊在外部，說明部分在內(nèi)部，應(yīng)該保留3個(gè)邊界點(diǎn)形成的三角形；如果網(wǎng)格3條邊或4條邊都在內(nèi)部，說明網(wǎng)格是內(nèi)部網(wǎng)格。

第3種是4個(gè)頂點(diǎn)都在邊界上的情況。如果網(wǎng)格中心點(diǎn)在內(nèi)部，說明是內(nèi)部網(wǎng)格。

2.5 需要變成三角形的網(wǎng)格

需變?yōu)槿切蔚木W(wǎng)格分為如下幾種。

第1種是上文2.4的第2種情況中提到的，此時(shí)應(yīng)將內(nèi)部點(diǎn)去掉，留下其余3個(gè)頂點(diǎn)變成三角形網(wǎng)格。

第2種是有且只有1個(gè)頂點(diǎn)在外部，即其它3個(gè)頂點(diǎn)可能在內(nèi)部或者邊界上。此時(shí)要做判斷，如果不是凹角附近的網(wǎng)格，則將這個(gè)外部點(diǎn)以及與它相鄰的2個(gè)點(diǎn)構(gòu)成的三角形，作為此處的網(wǎng)格，如圖6（a）所示。

第3種是2個(gè)對角的頂點(diǎn)在外，另外2個(gè)對角的頂點(diǎn)在內(nèi)。此時(shí)需要分成2個(gè)三角形。

后兩種情況變成三角形網(wǎng)格之后，需要對網(wǎng)格連接的合法性重新檢查，避免三角形網(wǎng)格穿過外形的內(nèi)部，如圖6（b）所示，頂點(diǎn)按照箭頭方向調(diào)整過后，避免了右上的三角形穿過外形內(nèi)部。

3 算法的分析

3.1 網(wǎng)格質(zhì)量

在光滑邊界的區(qū)域上，根據(jù)主要控制點(diǎn)的性質(zhì)，2個(gè)主要控制點(diǎn)的距離最近為[22]倍背景網(wǎng)格邊長。在尖點(diǎn)附近，如果四邊形網(wǎng)格的4個(gè)頂點(diǎn)中包含有尖點(diǎn)，雖無法保證最短邊的長度，但一個(gè)網(wǎng)格單元中的尖點(diǎn)只有一個(gè)，所以其它3條邊的長度是大于[22]倍邊長的。

如果三角形網(wǎng)格有一個(gè)頂點(diǎn)就是尖點(diǎn)，最小邊的長度至少為[12]倍邊長。所以這一方法在一定程度上保證了網(wǎng)格的質(zhì)量。

3.2 網(wǎng)格數(shù)量

本文的方法每個(gè)正方形網(wǎng)格最多分成2個(gè)三角形網(wǎng)格。網(wǎng)格生成的總數(shù)量不超過原始網(wǎng)格的2倍。且本文的方法僅僅在外形的邊界上至邊界外2層的范圍內(nèi)進(jìn)行計(jì)算，將需要計(jì)算的網(wǎng)格數(shù)量減少到最小。

3.3 時(shí)間復(fù)雜度分析

設(shè)邊界附近的網(wǎng)格數(shù)量為n。

2.1至2.5中的所有操作都是依賴于網(wǎng)格數(shù)量的，并且都是O（n）的復(fù)雜度。但因?yàn)?.3～2.5的步驟中都需要進(jìn)行分情況判斷，可能復(fù)雜度的常系數(shù)比較高。

對比Delaunay三角剖分的方法，因?yàn)榱鲌龅膬?nèi)外邊界有約束，所以這里需要用到約束的Delaunay三角剖分方法，時(shí)間復(fù)雜度最少為O（nlogn）。因此，只要n足夠大，本文的方法就會(huì)比帶約束的Delaunay三角剖分的方法快。

4 算法實(shí)例

5 結(jié)論

在二維復(fù)雜外形的邊界生成網(wǎng)格時(shí)，本文的方法根據(jù)外形與網(wǎng)格邊界的交點(diǎn)，分情況討論，制定出一套網(wǎng)格調(diào)整的方案，時(shí)間復(fù)雜度低于通常使用的Delaunay三角剖分方法。雖然網(wǎng)格的質(zhì)量達(dá)不到最優(yōu)的，但因?yàn)榻稽c(diǎn)的選取策略保證了網(wǎng)格點(diǎn)之間的距離，所以網(wǎng)格的質(zhì)量通常不會(huì)很差。

本文方法的優(yōu)勢是在網(wǎng)格質(zhì)量可以接受的限度內(nèi)，用復(fù)雜度O（n）的算法快速生成網(wǎng)格。

該方法為復(fù)雜外形的網(wǎng)格生成帶來了新的思路，今后還可以向在自適應(yīng)加密以及三維網(wǎng)格的方向推廣。

參考文獻(xiàn)：

[1] 齊學(xué)義，楊帆，齊沖.貼體坐標(biāo)網(wǎng)格生成技術(shù)的研究[J].工程熱物理學(xué)報(bào). 2001， 22（S1）：29-32.

[2] L.Paul Chew. Constrained delaunay triangulations [J]. Algorithmica. 1989， 4（1-4）：97-108.

[3] 楊名生，王冬.有限元網(wǎng)格全自動(dòng)生成的四分法[J].計(jì)算結(jié)構(gòu)力學(xué)及其應(yīng)用， 1989， 06（04）：61-68.

[4] 張來平，張涵信，高樹椿.矩形/非結(jié)構(gòu)混合網(wǎng)格技術(shù)及在二維/三維復(fù)雜無粘流場數(shù)值模擬中的應(yīng)用[J]. 空氣動(dòng)力學(xué)學(xué)報(bào)， 1998， 16（1）：79-88.

[5] W.-J. Ying， W.-C. Wang， A kernel-free boundary integral method for implicitly defined surfaces， Journal of Computational Physics 252 （2013） 606–624.

[6] 楊欽，劉瑞剛，孟憲海，等. 二維復(fù)雜限定Delaunay三角化算法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)， 2007， 19（02）：145-150.