自動生成四邊形網(wǎng)格的方法及其在數(shù)值模擬中的應(yīng)用

修榮榮,徐明海,黃善波

(1.中國石油大學(xué)期刊社,山東東營257061;2.中國石油大學(xué)儲運(yùn)與建筑工程學(xué)院,山東青島266555)

自動生成四邊形網(wǎng)格的方法及其在數(shù)值模擬中的應(yīng)用

修榮榮1,2,徐明海2,黃善波2

(1.中國石油大學(xué)期刊社,山東東營257061;2.中國石油大學(xué)儲運(yùn)與建筑工程學(xué)院,山東青島266555)

采用間接方法生成四邊形網(wǎng)格,首先利用改進(jìn)的兩點(diǎn)前沿推進(jìn)法把計(jì)算區(qū)域剖分成三角形網(wǎng)格,然后采用插點(diǎn)和細(xì)分的技術(shù)生成單元全部是四邊形的網(wǎng)格,通過邊互換、刪點(diǎn)和局部插點(diǎn)技術(shù)進(jìn)一步光滑平順,得到適用于數(shù)值計(jì)算的網(wǎng)格。剖分結(jié)果表明,該方法能夠在任意二維平面區(qū)域內(nèi)自動生成全四邊形網(wǎng)格,并能生成光滑過渡的局部加密網(wǎng)格和貼體性較好的邊界層網(wǎng)格。該方法具有算法簡單,計(jì)算量少的特點(diǎn)。利用所生成的網(wǎng)格對計(jì)算傳熱學(xué)中的典型算例-方腔自然對流進(jìn)行求解,計(jì)算結(jié)果與基準(zhǔn)解吻合,網(wǎng)格質(zhì)量能夠滿足數(shù)值分析計(jì)算的要求。

二維區(qū)域;全四邊形網(wǎng)格;三角形合并;變密度網(wǎng)格;數(shù)值模擬

在二維區(qū)域生成非結(jié)構(gòu)化網(wǎng)格的方法中,三角形網(wǎng)格的自動生成算法最為成熟,但有限元計(jì)算表明,在相同網(wǎng)格步長下,三角形單元的計(jì)算精度不如四邊形單元[1]。Aziz等[1]研究表明,有限容積法數(shù)值計(jì)算結(jié)果的精度與單元形狀關(guān)系不大,與網(wǎng)格步長關(guān)系則極為顯著。在同樣網(wǎng)格節(jié)點(diǎn)分布的情況下,三角形單元的數(shù)目是四邊形單元的兩倍,因此,最終形成的線性代數(shù)方程組的階數(shù)也是兩倍。這樣,無論是存儲還是求解,都是不經(jīng)濟(jì)的。因此,為了在一定計(jì)算工作量的條件下盡量提高精度,提出了采用非結(jié)構(gòu)化四邊形網(wǎng)格計(jì)算思想,進(jìn)而導(dǎo)致了四邊形網(wǎng)格生成的技術(shù)問題。生成四邊形網(wǎng)格的方法很多[2-8],按照單元生成方法大致分為兩類:直接法和間接法。直接法[2-3,8]包括幾何分割法和前沿推進(jìn)法,其生成的網(wǎng)格質(zhì)量高,不規(guī)則節(jié)點(diǎn)少,但在網(wǎng)格生成的過程中需要進(jìn)行交叉檢驗(yàn)和局部光滑(如paving方法[3]),網(wǎng)格生成速度較慢。相對而言,間接法的操作都是局部的,且不必進(jìn)行交叉檢測,因而網(wǎng)格生成的速度快。間接法[4-5,7]包括插入節(jié)點(diǎn)法和合并三角形法[5]。前者能把三角形網(wǎng)格完全轉(zhuǎn)化為四邊形網(wǎng)格,但插入的節(jié)點(diǎn)導(dǎo)致大量不規(guī)則單元,網(wǎng)格單元質(zhì)量較差;后者提高了單元的質(zhì)量,但不能保證把所有的三角形單元轉(zhuǎn)換成四邊形單元。為克服殘留三角形問題,文獻(xiàn)[6]對其作了進(jìn)一步的改進(jìn),可以把三角形網(wǎng)格完全轉(zhuǎn)化為四邊形網(wǎng)格,且可實(shí)現(xiàn)局部加密,保證網(wǎng)格質(zhì)量,但由于用布點(diǎn)法生成三角形網(wǎng)格,生成的變密度網(wǎng)格不夠光滑。筆者基于兩點(diǎn)前沿推進(jìn)法[9]首先將計(jì)算區(qū)域剖分成三角形網(wǎng)格,并通過函數(shù)直接給出背景網(wǎng)格的信息來控制網(wǎng)格的疏密,然后再將三角形網(wǎng)格合并成四邊形網(wǎng)格,以改進(jìn)網(wǎng)格質(zhì)量。

1 網(wǎng)格生成算法

生成四邊形網(wǎng)格的算法基于三角形網(wǎng)格生成過程,經(jīng)過拆分和局部優(yōu)化,把三角形網(wǎng)格轉(zhuǎn)化為四邊形網(wǎng)格。

1.1 三角形網(wǎng)格的生成

采用改進(jìn)的兩點(diǎn)前沿推進(jìn)算法生成三角形網(wǎng)格[9]。設(shè)邊界上指定點(diǎn)處的四邊形網(wǎng)格劃分尺寸為h,則在三角形網(wǎng)格劃分時(shí)離散尺寸為2h[6]。為克服文獻(xiàn)[6]由于采用布點(diǎn)法控制網(wǎng)格步長而導(dǎo)致變密度網(wǎng)格不夠光滑的缺點(diǎn),引入了通過函數(shù)直接給出背景信息的處理方法[10],先生成變密度的三角形網(wǎng)格,從而得到光滑過渡的變密度網(wǎng)格,使人為的工作量和計(jì)算量大大減少。方法的基本原理是利用電荷產(chǎn)生的勢函數(shù)隨距離衰減的規(guī)律作為網(wǎng)格步長變化的控制函數(shù),再結(jié)合電荷產(chǎn)生的勢場可疊加的原理,通過多個(gè)點(diǎn)的網(wǎng)格背景信息,實(shí)現(xiàn)網(wǎng)格步長的任意加密,通過適當(dāng)布置一定數(shù)量和種類的基本解作為網(wǎng)格步長控制函數(shù),替代傳統(tǒng)的背景網(wǎng)格,可以節(jié)省網(wǎng)格步長控制搜索時(shí)間,進(jìn)而提高網(wǎng)格生成速度。

1.2 四邊形網(wǎng)格的生成

四邊形網(wǎng)格生成的步驟如下:

(1)合并三角形。對三角形單元的集合T中任一單元ABC,與其共邊的三角形最多有3個(gè),分別計(jì)算三角形ABC與它們合并后所得四邊形的網(wǎng)格質(zhì)量(記為β),取最大值所對應(yīng)的三角形為候選合并三角形。若β值大于等于0,則執(zhí)行合并操作,將四邊形加到四邊形集合Q中;否則不予合并,將三角形ABC加到殘留三角形集合Tr中。重復(fù)此過程,直至集合T為空集。

(2)三角形的細(xì)分。由于三角形合并成四邊形時(shí)總有殘留三角形存在,這些三角形無法與其相鄰三角形合并成四邊形。為生成四邊形,在三角形的形心處插入一點(diǎn),與其三邊中點(diǎn)相連,形成3個(gè)四邊形,但這樣破壞了其相鄰四邊形的拓?fù)浣Y(jié)構(gòu),為此需要將三角形合并而成的四邊形也要細(xì)分一次,以滿足網(wǎng)格的拓?fù)湟蟆?/p>

(3)四邊形的細(xì)分。對集合Q中的每一個(gè)四邊形ABCD,將其形心O與各邊中點(diǎn)連接,形成四個(gè)小四邊形。

對于剩余的三角形處理過程如圖1所示?？梢钥闯?這種方法可以把三角形網(wǎng)格完全轉(zhuǎn)化為四邊形網(wǎng)格。因?yàn)樵瓉聿介L為2h,所以細(xì)分后步長恰好滿足要求。

圖1 剩余三角形的細(xì)分及網(wǎng)格處理過程Fig.1 Subdivision of residual triangles and mesh proceesing procedure

2 網(wǎng)格質(zhì)量的改進(jìn)

2.1 四邊形的質(zhì)量判斷準(zhǔn)則

如圖2(a)所示,按逆時(shí)針走向的四邊形ABCD被對角線AC分成兩個(gè)三角形ABC和ACD,其質(zhì)量用α表示,分別為同時(shí),四邊形又可被對角線BD分成兩個(gè)三角形ABD和BCD,其α值分別為按降序重新排列成為α1,α2,α3,α4,即α1≥α2≥α3≥α4,則四邊形ABCD的質(zhì)量判斷準(zhǔn)則[11]為

其中,按逆時(shí)針走向的三角形ABC(圖2(b))的質(zhì)量按下式計(jì)算:

式中,n為三角形單位法向。α的最大值為1,所對應(yīng)的三角形為等邊三角形,α值越小,三角形質(zhì)量越差。

圖2 凸四邊形、三角形和凹四邊形示意圖Fig.2 Convex quadrilateral,triangle and concave quadrilateral

對于凸四邊形(圖2(a)),β的值在0和1之間,矩形的值最大為1。凹四邊形(圖2(c))的β值小于0。因此,可認(rèn)為β的值越小四邊形質(zhì)量越差。

2.2 網(wǎng)格的優(yōu)化光順

最終網(wǎng)格質(zhì)量的優(yōu)劣對數(shù)值結(jié)果的精度和收斂性至關(guān)重要,生成高質(zhì)量網(wǎng)格的一個(gè)重要措施就是進(jìn)行網(wǎng)格后處理,即網(wǎng)格優(yōu)化[12]。網(wǎng)格的優(yōu)化分兩步進(jìn)行:局部優(yōu)化[13]和整體光順。

(1)局部優(yōu)化。由于三角形內(nèi)插入一點(diǎn)形成的四邊形的單元質(zhì)量較差,故在經(jīng)過網(wǎng)格光順處理后,還需要進(jìn)一步處理以提高單元質(zhì)量。四邊形網(wǎng)格中的每個(gè)內(nèi)節(jié)點(diǎn)的最佳相鄰單元數(shù)為4,這樣與其相關(guān)的單元接近于矩形。如果相鄰單元數(shù)(設(shè)為Ne)遠(yuǎn)遠(yuǎn)大于或小于4,即那么以該內(nèi)節(jié)點(diǎn)為頂點(diǎn)的四邊形必然有不規(guī)則的。為此,對滿足下列條件的單元進(jìn)行刪除:一個(gè)單元的兩相對節(jié)點(diǎn)A、C都是可動節(jié)點(diǎn),且相鄰單元數(shù)都為3,如圖3所示。

圖3 單元?jiǎng)h除Fig.3 Element deletion

(2)網(wǎng)格的整體光順。采用Laplace法進(jìn)行光順處理可以改善網(wǎng)格的質(zhì)量。光順的本質(zhì)是移動內(nèi)部網(wǎng)格節(jié)點(diǎn)的位置,使其處于有共享此點(diǎn)的四邊形單元所組成的多邊形的中心。光順處理是以改善局部最差四邊形為標(biāo)準(zhǔn),即對任一內(nèi)部節(jié)點(diǎn)D(x,y),與其相連的四邊形最小β值為βmin,假設(shè)經(jīng)光順后的節(jié)點(diǎn)為D＊(x＊,y＊),與其相連的四邊形最小β值為則用D＊(x＊,y＊)代替D(x,y),否則不做處理。

3 剖分算例

對單連通區(qū)域和多連通區(qū)域等不規(guī)則區(qū)域進(jìn)行了剖分,經(jīng)過光滑后的剖分結(jié)果如圖4所示。從圖中可以看出,可以生成均勻網(wǎng)格、光滑過渡的局部加密網(wǎng)格和貼體性較好的邊界層網(wǎng)格。

圖4 剖分實(shí)例Fig.4 M esh generation results

4 數(shù)值模擬算例

利用上述方法生成的網(wǎng)格,計(jì)算了已有比較一致結(jié)果的方腔自然對流及半圓形空腔內(nèi)的自然對流問題,對剖分的網(wǎng)格進(jìn)行了考核,并與三角形網(wǎng)格進(jìn)行了對比。

4.1 方腔自然對流

采用網(wǎng)格步長相同的三角形網(wǎng)格與四邊形網(wǎng)格分別對Ra=103,104,105,106進(jìn)行求解。三角形網(wǎng)格的節(jié)點(diǎn)數(shù)為1598,單元數(shù)為3050;四邊形網(wǎng)格的節(jié)點(diǎn)數(shù)為1977,單元數(shù)為1896。

4.1.1 計(jì)算精度的對比

將兩種網(wǎng)格計(jì)算所得的平均努塞爾(Nusselt)數(shù)及其對應(yīng)的位置與文獻(xiàn)[14]認(rèn)為最好的解進(jìn)行比較,比較結(jié)果見表1。其中,局部努塞爾數(shù)(N u)與熱壁上的平均努塞爾數(shù)按下式進(jìn)行計(jì)算:

表1中數(shù)據(jù)表明,在單元邊長基本相同的情況下,四邊形網(wǎng)格的計(jì)算精度優(yōu)于三角形網(wǎng)格,尤其是最大、最小努塞爾數(shù)和對應(yīng)的位置方面。

表1 努塞爾數(shù)計(jì)算結(jié)果與基準(zhǔn)解的比較Table 1 Comparison of calculated results of Nusselt number and reference results

4.1.2 收斂特性的對比

收斂特性曲線如圖5所示。圖5(a)中三角形網(wǎng)格的節(jié)點(diǎn)數(shù)為710,單元數(shù)為1322,四邊形網(wǎng)格的節(jié)點(diǎn)數(shù)為753,單元數(shù)為704;圖5(b)中三角形網(wǎng)格的節(jié)點(diǎn)數(shù)為1598,單元數(shù)為3 050,四邊形網(wǎng)格的節(jié)點(diǎn)數(shù)為1637,單元數(shù)為1 564。由圖中可以看出,網(wǎng)格步長相同時(shí),四邊形網(wǎng)格的收斂速率明顯高于三角形網(wǎng)格。

圖5 收斂速率的比較Fig.5 Comparison of convergence rate

4.1.3 所耗機(jī)時(shí)的對比

兩種網(wǎng)格的網(wǎng)格節(jié)點(diǎn)數(shù)與占用機(jī)時(shí)的關(guān)系曲線見圖6。從圖中可以看出,網(wǎng)格節(jié)點(diǎn)數(shù)相同時(shí),四邊形網(wǎng)格所耗機(jī)時(shí)遠(yuǎn)小于三角形網(wǎng)格,而且網(wǎng)格節(jié)點(diǎn)數(shù)越多,節(jié)省時(shí)間越明顯。

4.2 半圓形空腔內(nèi)的自然對流

采用四邊形網(wǎng)格(節(jié)點(diǎn)數(shù)為811,單元數(shù)為758)計(jì)算出Ra=104,105時(shí)熱邊界上努塞爾數(shù)的分布曲線(圖7),圖中熱壁上X值的范圍為0.15≤X≤0.85。可以看出,本文的計(jì)算結(jié)果與文獻(xiàn)[15]是一致的。

5 結(jié)束語

本文中提出的四邊形網(wǎng)格生成方法具有計(jì)算量少、程序量少且易編制的特點(diǎn),能夠在任意二維平面區(qū)域生成全四邊形網(wǎng)格,可以靈活地實(shí)現(xiàn)網(wǎng)格的局部加密,特別適合自適應(yīng)網(wǎng)格計(jì)算的要求。利用所生成的四邊形網(wǎng)格對計(jì)算傳熱學(xué)中的兩個(gè)典型算例的求解結(jié)果與基準(zhǔn)解一致,說明網(wǎng)格質(zhì)量能夠適應(yīng)數(shù)值分析計(jì)算的要求。

[1] AZIZ K.Reservoir simulation grids:opportunities and problems[R].SPE 25233,1993.

[2] 趙熠,趙建軍,張新訪.前沿法生成四邊形網(wǎng)格的改進(jìn)方法[J].計(jì)算機(jī)工程與應(yīng)用,2002(9):64-66,98.ZHAO Yu,ZHAO Jian-jun,ZHANG Xin-fang.Generating quadrilateral mesh based on the improved advancing front method[J].Computer Engineering and Applications,2002(9):64-66,98.

[3] BLACKER TD,MARK S S.Paving:a new approach to automated quadrilateral mesh generation[J].International Journal For Numerical Methods in Engineering,1991,32(4):811-847.

[4] 馮道雨,陳尚法,陳勝宏.復(fù)雜區(qū)域生成四邊形網(wǎng)格的一種改進(jìn)方法[J].巖土力學(xué),2004,25(6):917-921.FENG Dao-yu,CHEN Shang-fa,CHEN Sheng-hong.A new quadrilateral mesh generation method based on advancing front technique[J].Rock and Soil Mechanics,2004,25(6):917-921.

[5] LO S H.Generating quadrilateral elements on plane and over curved surfaces[J].Comput Struct,1989,31:421-426.

[6] 顧元憲,馬正陽,關(guān)振群.平面任意區(qū)域四邊形網(wǎng)格自動生成的一種方法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),1998,10(5):432-439.GU Yuan-xian,MA Zheng-yang,GUAN Zhen-qun.A method of automatic generation of quadrilateral element meshes on arbitrary plane domain[J].Journal of Computer Aided Design&Computer Graphics,1998,10(5):432-439.

[7] 劉晶,聶玉峰,蘇少普.四邊形網(wǎng)格間接生成方法[J].計(jì)算機(jī)工程與應(yīng)用,2010,46(2):44-47.L IU Jing,N IE Yu-feng,SU Shao-pu.Indirect method of quadrilateral mesh generation[J].Computer Engineering and Applications,2010,46(2):44-47.

[8] 李毅,鮑勁松,金燁,等.二維域多約束四邊形有限元網(wǎng)格生成算法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2008,20(4):488-493.L I Yi,BAO Jin-song,J IN Ye,et al.Quadrilateralmeshgeneration algorithm for planar domain with multi-constraints[J].Journal of Computer Aided Design&Computer Graphics,2008,20(4):488-493.

[9] 修榮榮,徐明海,黃善波,等.二維區(qū)域三角形的一種改進(jìn)的前沿推進(jìn)法[J].石油大學(xué)學(xué)報(bào):自然科學(xué)版,2003,27(5):73-75,80.X IU Rong-rong,XU Ming-hai,HUANG Shan-bo,et al.An improved triangulation method for arbitrary two-dimensional planar domain using advancing method[J].Journal of the University of Petroleum,China(Edition of Natural Science),2003,27(5):73-75,80.

[10] 武潔,馮晉利.三角形網(wǎng)格生成方法中一種提供背景信息的方法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),1999,11(4):300-303.WU Jie,FENG Jin-li.The background inforrnation in advancing-front method for the triangle meshing[J].Journal of Computer Aided Design&Computer Graphics,1999,11(4):300-303.

[11] LEE C K,LO S H.A new scheme for the generation of a graded quadrilateral mesh[J].Comput Struct,1994,52:847-857.

[12] 陳立崗,鄭耀,陳建軍.全四邊形有限元網(wǎng)格的拓?fù)鋬?yōu)化策略[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2007,19(1):78-83.CHEN Li-gang,ZHENG Yao,CHEN Jian-jun.Topological improvement for quadrilateral finite element meshes[J].Journal of Computer Aided Design&Computer Graphics,2007,19(1):78-83.

[13] ZHU J Z,ZIENKIEW ICZ O C,H INTON E,et al.A new approach to the development of automatic quadrilateral mesh generation[J].Int J Numer Meth Engng,1991,29:1551-1567.

[14] Devahl DAV IS G,JONES I P.Natural convection in a square cavity,A comparison Exercise[M]//LEW IS R W,MORGAN K,SCHREFLER B A.Numerical methods in thermal problems.Swansea,U K:Pineridge Press,1981:552-572.

[15] VO ILKER S,BURTON T,VANKA S P.Finite-volume multigrid calculation of natural-convection flows on unstructured grids[J].Numerical Heat Transfer(PartB),1996,30:1-22.

(編輯 韓國良)

Automatic generation method of quadrilateral meshes and its application to numerical simulation

XIU Rong-rong1,2,XU Ming-hai2,HUANG Shan-bo2

(1.Periodical Office of China University of Petroleum,Dongying257061,China;2.College of Storage&Transportation and Architectural Engineering in China University of Petroleum,Qingdao266555,China)

Quadrilateral meshes on arbitrary 2D domain were generated by using indirect method.Firstly,triangle meshes were generated by improved two-point advancing method.Then by means of inserting nodes and subdivison technique,the mesh elements were composed with quadrilaterals completely.The edge interconversion,node elimation and inserting node locally were used to smooth meshes,and the meshes suitable to numerical calculation were obtained.The proposed method is characterized by its simplicity and ease of implementation.Mesh generation results demonstrate the robustness of the method.The method can also generate smooth graded quadrilateral mesh and body-fitted boundary layer meshes.The generated meshes were used to solve natural convection in a square cavity.The solution results agree well with the reference results.The quality of the generated meshes can be satisfied with the requirements of finite element numerical simulation.

two-dimensional domain;all-quadrilateral mesh;triangle merging;graded mesh;numerical simulation

TK 124;TP 391.41

A

10.3969/j.issn.1673-5005.2011.02.023

2010-04-20

國家自然科學(xué)基金項(xiàng)目(50904077);中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助項(xiàng)目(10CX05011A)

修榮榮(1977-),女(漢族),山東萊陽人,碩士,主要從事CDF/NHT技術(shù)研究。

1673-5005(2011)02-0131-06