參數曲面三角網格生成的改進波前法

王偉,樊宏周,席光

(西安交通大學能源與動力工程學院, 710049, 西安)

參數曲面三角網格生成的改進波前法

王偉,樊宏周,席光

(西安交通大學能源與動力工程學院, 710049, 西安)

為了消除基于波前法的有限元三角網格算法在參數曲面網格剖分過程中單元形狀映射畸變的問題,結合直接法和映射法各自的優點,提出了一種新的三角網格生成算法,即:對當前節點進行剖分,并在三維空間直接產生新節點且進行節點的合法性判斷,再將物理網格映射到參數空間形成參數域網格;對相鄰波前段形成的角度進行剖分,依據角度大小生成個數不等的單元,通過優先剖分銳角節點使波前段始終構成鈍角多邊形。經剖分算例表明:所提算法減少了節點合法性判斷內容和判斷次數,避免了重復剖分,取消了剖分結束算法,提高了網格剖分效率,生成了高質量的三角網格;僅需對網格排列情況的直觀分析,便可定性判斷三維曲面的空間曲率變化。該算法對葉片加工中振動分析、精密加工研究等具有指導意義。

有限元;曲面網格;波前法;三角剖分;映射畸變

曲面網格生成是有限元分析中幾何建模的關鍵,曲面網格可以直接作為膜結構、殼體等的有限元模型進行計算分析。三維實體的表面是由各種曲面和平面組成,三維實體表面離散的數據點及曲面網格質量將直接影響剖分算法的可靠性和生成體網格的質量[1]。所以,有限元、計算流體動力學和計算機圖形學分析中獲得高質量的曲面網格是必不可少的前提條件[2]。

曲面網格剖分算法中波前法(Advancing Front Technique,AFT)是曲面三角網格生成的重要方法,具有啟發性、局部適應性和單元可控性[3]。基于波前法的三角網格剖分算法涉及到映射法和直接法。映射法是運用波前法將參數空間的參數域進行剖分,并映射到三維曲面上生成曲面網格,其研究重點是克服曲面參數化不均勻導致的網格映射畸變。為了解決映射畸變問題,Cuillière先由曲面參數方程推導出描述映射畸變程度的解析式,再運用改進的波前法來剖分參數域,然后將參數域映射到曲面上以形成曲面網格[4-5]。這種由解析式對映射畸變的度量,其結果是準確的,但通用性略顯不足。Lee等在采用映射法剖分復雜三維曲面時,將整個曲面分割成若干小而簡單的曲面片,然后各自進行剖分,在參數域的剖分中運用黎曼度量剖分后得到各向異性網格,這樣在三維曲面上得到理想的三角網格[6-7]。可見,映射法在網格映射過程中的映射畸變是制約該方法生成高質量網格的關鍵。直接法是運用波前法直接在曲面上產生新節點并生成新單元,是波前法從二維空間向三維空間的擴展,其研究重點是在三維空間曲面上精確定位新節點,生成高質量單元,做高效的網格剖分。在三維空間生成新節點方面:Lo先在剖分區域構造輔助曲線,再采用線面求交的方法精確定位曲面上的新節點[8],該方法對剖分曲率變化平緩的空間曲面能夠得到高質量網格;Ito等在曲面的網格節點中采用插值的方法生成新節點,但新節點過度依賴曲面上網格節點的數量和分布情況[9]。可見,在新節點生成之前直接法缺少對區域有效的預判環節,從而導致新節點的產生具有盲目性,區域剖分出現重復進行,對新節點的合法性要進行大量判斷,包括相交判斷、包含判斷及距離判斷等,這種判斷耗用了大約80%的計算時間,既影響剖分效率,又影響網格質量。

本文運用映射法的思路在三維曲面上為每個節點建立正交標架,運用直接法的思路、采用投影算法直接在三維曲面上定位新節點。在節點的正交標架中,通過投影算法定位新節點,在三維空間中結合映射法、直接法直接對新節點進行合法性判斷。曲面剖分時先對相鄰波前段形成的角度進行剖分,再結合節點間的距離生成單元網格。為避免區域剖分的重復性,優先剖分銳角節點,從而減少了節點合法性判斷的次數,縮短了算法的運行時間,增強了算法的可靠性,提高了剖分效率和網格質量。

1 新節點生成

在生成三角網格的過程中,波前法將邊界上的網格點所形成的一系列有向線段的集合定義為波前,從波前出發,逐一與區域內部的點形成三角形單元,該三角形單元的新邊加入到波前的行列,生成該單元的邊則從波前行列中消除,如此不斷地向區域內部推進,直到波前的行列為空為止,此時生成的三角網格覆蓋整個區域,剖分結束。

在剖分網格的過程中,本文借鑒了映射法中的坐標系,即在三維曲面上為任意節點建立正交標架,并在2個空間中分別建立具有一一對應關系的坐標系,而在三維曲面上借鑒了直接法的節點生成方式直接定位新節點。所以,本文算法在2個空間共建立了3個坐標系,在三維曲面上的坐標系中運用投影算法直接產生新節點。這樣,在一個節點的產生過程中,本文算法集中了映射法和直接法各自的優點,巧妙避免了映射畸變的問題。

波前法的關鍵是新節點生成算法,其在待剖分的三維曲面上,根據不同部位的特點產生相應的新節點。對于遠離邊界約束的中心區域,應盡可能產生理想形狀的單元新節點;對于邊界附近的區域,受邊界、角點等因素的約束,應靈活生成不同形狀的單元新節點。

(1)

式中:ru、rv分別為曲面上點M沿u、v的單位切向量。

(a)參數空間 (b)物理空間

方程(1)表示線段MC在XL方向上的投影為0,在ZL方向上的投影為dMC,求解該方程可得相應參數的增量Δu和Δv。點c在參數空間的坐標為

(2)

在確定了當前剖分節點B之后,通過調整φ和dBC可以在三維曲面上任意定位新節點C,生成各種不同形狀的三角單元。理想單元顯然不能覆蓋角點、極點等特殊區域,這就需要通過合理調整φ和dBC靈活移動新節點的位置,以運用高質量的單元完全覆蓋所有區域。新節點生成算法利用了波前法的良好的局部適應性和靈活性來調整φ和dBC。

2 新單元生成

波前法生成三角網格是一種探索推進式的單元生成過程,其根據已有波前向區域內部延伸網格,在產生新節點的同時生成新單元,并隨時檢驗新節點的合法性。波前段的大小直接反映已生成單元的質量,同時影響后續將要生成的新單元的質量,對波前段大小進行嚴格的控制,可保證生成高質量的單元。當前剖分節點B的確定決定著網格剖分的順序,φ和dBC的確定決定著新生成單元的形狀。節點的合法性判斷占用了算法大部分的計算時間,所以壓縮判斷內容、減少判斷次數,才能提高剖分算法的效率。在新單元的生成過程中,將兩相鄰波前形成的角度定義為剖分角度θ,本文算法通過對θ的剖分,再配合線段長度,將生成合理的新單元。

2.1 波前段極值

三角單元的邊長會影響網格質量,波前段是連接已有單元和后續新單元的紐帶,其大小直接影響新單元的質量,本文設定了單元邊長的2個極值。

令單元邊長的極小值

(3)

單元邊長的極大值

(4)

式中:d為理想單元邊長。

將網格單元的邊長嚴格控制在2個極值下,不會生成單個邊長過小的銳角單元和單個邊長過大的鈍角單元。

2.2 三角網格生成

區域外邊界初始化的波前以逆時針方向為正向,區域內邊界初始化的波前以順時針方向為正向,這樣可確保剖分區域始終位于波前正方向的左側。每個波前節點都有相應的剖分角度θ,隨著波前的不斷更新,θ需要適時更新。以節點P3(見圖2)為當前剖分節點,根據θ的大小判斷局部區域的形態,結合節點間的長度生成新節點,產生新單元。θ可以分為以下幾種。

(a)剖分角度過小 (b)線段P2P4長度適當

(c)線段P2P4長度過短 (d)線段P2P4長度過長

(e)一個新節點 (f)第1個新節點

(1)當dmin≤‖P2P4‖≤dmax時,如圖2b所示,表明該區域適于生成一個新單元P2P3P4,此時φ=θ,應在波前節點列中取消節點P3;

(2)當‖P2P4‖

(5)

如果‖d36-d‖>0.1d,則令d36=d,以保證新生成的波前段盡量接近d,進而保證后續新單元的大小盡量接近d。在此,算法對線段的長度范圍進行了嚴格限制,使波前段的長度盡量接近d,以確保優化后續的剖分形態。由新節點P6生成2個新單元P2P3P6和P3P4P6后,波前節點列中以節點P6替代節點P3。在三角剖分過程中,剖分角度為鈍角時,表明該節點區域開闊,可以生成給定要求的新單元,且不易發生重復剖分,也不存在優先剖分,可按照算法要求的剖分順序確定當前的剖分節點。

(6)

根據單元數將θ均分,由此得到φ,順時針依次解三角形,運用新節點的生成算法依次生成相應的新節點和新單元,直至將θ剖分完畢為止。

本文算法創新性地角度剖分,為波前形態的描述提供了確切的依據,解決了傳統波前法對于圖3所示的部分波前P1P2P3及新節點P4、P5中,相鄰單元容易相交(見圖3a)或距離過近(見圖3b),從而造成剖分失敗、算法復雜度增加的問題。

(a)相鄰單元相交 (b)相鄰單元距離過近

從上述剖分角度知,本文算法不必對區域角點單獨進行分析,區域角點只是普通的剖分角度,可見按角度大小對波前段進行分類并分別剖分,可減少運算時間,降低算法的復雜度,提高網格的質量。

2.3 節點合法性判斷

經過上述判斷,新節點P6可能會與自身波前在某一個節點處連接起來,這樣一個波前被分割成了2個波前;新節點P6也可能會與其他波前在某一個節點處連接起來,這樣2個波前被合并成一個波前。總之,隨著三角剖分的推進,波前數越來越少,波前節點數越來越少,當波前只有一個且其中的節點數為0時,剖分結束,三角網格覆蓋整個區域。

從上述節點合法性判斷知,本文算法在節點的合法性判斷方面僅需考慮新節點與波前的其他節點間的距離、新單元與其他波前段是否相交的問題,相比傳統的AFT算法,減少了判斷的內容和次數,縮短了計算時間,且可一次性高質量完成剖分任務。

2.4 算法流程

從上述算法流程知,本文算法面對的剖分對象是節點,剖分可以進行到波前段為空,即波前節點為2個或3個時,剖分結束,直接結束程序卻不需要進行繁瑣的剖分結束判斷。

3 網格質量系數

文獻[10]對三角網格質量的衡量準則進行了詳細分析,本文分別采用形狀質量系數和整體質量系數對剖分的三角網格從2個不同的角度進行了評估。對于單元ABC的形狀質量系數

(7)

區域網格的形狀質量系數

(8)

(9)

式中:NT為網格總數;α1、α2為相同衡量準則下對網格做出共同評價的一對相關系數。

依據文獻[6]可以得到網格單元的整體質量系數μ及區域網格的整體質量系數μ1和μ2。μ1、μ2作為相同衡量準則下的一對相關的系數共同對網格質量做出評價。

形狀質量系數是評價網格單元形狀的系數,等邊三角單元的形狀質量系數為1。由于本文算法在三維空間直接產生了新節點和新單元,不會發生映射畸變,所以運用該系數可以準確衡量算法生成單元的質量。

從網格質量衡量準則知,本文算法針對區域的不同部分,可以靈活生成不同形狀的三角單元,所以運用所提系數可以準確衡量覆蓋整個區域網格的整體質量,避免了文獻[6]只有在保證單元的大小與理想單元接近時才能保證區域網格整體質量的問題。

4 剖分算例

運用本文算法剖分了平面上的正方形,比較了2種剖分方式的優劣;剖分了平面上的多聯通區域,以檢驗算法對多個波前剖分的能力;通過平面區域剖分驗證了算法對角度剖分的可行性及生成網格的質量。在三維空間進行了圓柱面剖分,所得三角化網格與解析計算結果完全一致,如圖4所示。

圖4 圓柱曲面網格

4.1 二維空間方形剖分

圖5a運用波前法以方形的四條邊界為波前進行了環形推進剖分,在四條邊界附近生成了高質量的網格,在區域內部波前交界處形成如水流交匯的網格,最終在方形中心結束剖分。圖5b運用波前法以上下2個相對的邊界為波前進行剖分,在初始的2個邊界附近和區域內部生成布局整齊的高質量網格,在另外2個邊界附近和剖分結束區域形成明顯的波前剖分相交的痕跡,生成排列復雜的網格。由此可知,波前法可以根據剖分區域的要求靈活選擇剖分方式。如果剖分區域的邊界是后續計算分析的關鍵部位,則可采取所有邊界為波前進行環形剖分;如果剖分區域的中間部分為后續分析的重點,則可采取兩相對較長邊界為波前的剖分方式。總之,初始波前附近的網格,由于是剖分的起始部分,所以網格質量一般都比較高,排列較為整齊。

(a)按邊界環形推進 (b)按對邊推進

4.2 二維多聯通區域剖分

如圖6所示,運用波前法剖分了平面多聯通區域,共生成了5363個三角網格單元,單元的質量系數α1=0.988 3、α2=0.987 9、μ1=0.9656、μ2=0.956 2(見圖6a)。該剖分以5個邊界分別為初始波前同時推進剖分,而在5個邊界附近明顯存在按層推進剖分的高質量網格。隨著剖分的推進,5個波前相互合并,形成類似水流交匯的網格排列,最終合并為一個波前,直至剖分結束。網格的形狀質量系數較高,表明所得網格的形狀較好;整體質量系數偏低,表明該層網格中有過大、過小的網格,該類網格恰恰出現在波前交匯區域。由于區域中有5個初始波前最終合并為1個波前,經4次波前合并后有4個明顯的波前交界線,其會導致網格的整體質量系數偏低。由此可見,本文算法可高質量地剖分多聯通區域,所得網格的排列情況也可客觀地反映多聯通區域的構成特點。

(a)平面聯通區域三角網格 (b)正螺旋面三角網格

4.3 正螺旋面的三角剖分

在三維空間運用波前法剖分了正螺旋面,共生成5500個三角網格單元,網格質量系數α1=0.9858、α2=0.9854、μ1=0.963 9、μ2=0.961 3(見圖6b)。三維空間中該曲面扭曲,空間曲率變化大。三角網格剖分時首先在曲面上建立曲紋坐標,以曲面的4個邊界為初始波前,按層推進剖分。在曲面的邊界附近存在具有明顯層次排列的網格,網格剖分至空間曲率變化大的曲面中心時,三角網格的排列失去層次,單元分布較為復雜,但由質量系數可知,所得三角網格依然是高質量的,表明本文算法可以對空間扭曲的三維曲面進行高質量的三角剖分。

4.4 三維自由曲面葉片的剖分

圖7為離心葉輪三維葉片的剖分結果。離心葉輪葉片的空間形狀為三維自由扭曲曲面,其網格包含了3 881個三角單元,由葉片的4個邊界附近的區域得到的排列整齊的三角網格保留了按層推進剖分的痕跡。整個葉片的網格布局層次分明,單元排列整齊。葉片的網格質量系數α1=0.992 0、α2=0.991 8、μ1=0.980 6、μ2=0.979 2。由圖7c可知,2條質量系數曲線幾乎重合,表明2個系數衡量的網格質量基本一致。形狀質量系數顯示出單個單元的形狀較好,可見采用本文算法在每個新節點的生成過程中都能精確定位新節點。

(a)四邊形結構化網格

(b)三角非結構化網格

(c)網格質量系數

綜上所述,本文算法在剖分三維扭曲曲面時所得的三角網格在整體布局上層次分明,單元排列整齊,網格整體質量較高,凸顯了直接法和映射法的優點同時避免了各自的不足。

通過本文算法剖分所得網格還可以定性判斷三維曲面上的曲率變化情況。如圖7b所示,在葉片的出口部分,網格層次分明、清晰,可以清楚地分辨出波前交匯處的網格,結合葉輪葉片的形狀可知,這里正是葉片曲率變化微小的區域。葉片的進口中心區域(圖7b中標出的區域),網格排列復雜,層次不甚分明,這里正是葉片曲率變化劇烈、空間扭曲大的部位。可見,網格排列整齊與否和曲面的曲率變化存在著明顯的對應聯系。圖7b中標出的區域位于葉片進口靠近葉根側,表明該葉片在進口中心靠近葉根區域的曲率變化更大,靠近葉頂區域的曲率變化稍顯平緩。由上述分析可知,僅僅通過對三維曲面剖分所得網格的直觀分析,可定性判斷曲面的空間曲率變化,這對后續的葉片加工中的振動分析、葉片精密加工研究等環節具有指導作用。

5 結 論

本文對波前法三角網格剖分進行了改進,結合直接法的節點生成方式和映射法的坐標體系,提出了一種新的三角網格生成算法。該算法既解決了映射法中的映射畸變問題,又解決了直接法的剖分盲目性問題。在節點的正交標架中,運用投影算法產生的新節點可以生成任意形狀的三角單元,充分發揮了波前法的局部靈活性的優勢。通過對角度的剖分,推進了網格的生成,通過對銳角節點的優先剖分及局部區域形態的判斷,避免了直接法剖分的盲目性,減少了節點合法性判斷的次數,提高了剖分效率。

采用本文算法僅通過對網格排列情況的直觀分析,即可定性判斷三維曲面的空間曲率變化情況。剖分算例表明,本文算法可靠、高效,所生成的網格層次分明、排列整齊,是一種符合有限元分析要求的高質量網格生成算法。

[1] 黃曉東, 丁問司, 杜群貴.基于波前法的參數曲面有限元網格生成算法 [J].計算機輔助設計與圖形學學報, 2010, 22(1): 51-59.HUANG Xiaodong, DING Wensi, DU Qungui.Parametric surface mesh generation based on advancing front technique [J].Journal of Computer Aided Design & Computer Graphics, 2010, 22(1): 51-59.

[2] 樊文剛, 李建勇, 黃澤華, 等.多點切觸加工在復雜凸曲面中的應用 [J].西安交通大學學報, 2012.46(3): 53-58.FAN Wengang, LI Jianyong, HUANG Zehua, et al.Application of multi-point contact machining to convex sculptured surface [J].Journal of Xi’an Jiaotong University, 2012, 46(3): 53-58.

[3] 關振群, 宋超, 顧元憲, 等.有限元網格生成方法研究的新進展 [J].計算機輔助設計與圖形學學報, 2003, 15(1): 1-13.GUAN Zhenqun, SONG Chao, GU Yuanxian, et al.Recent advancing of research on finite element mesh generation methods [J].Journal of Computer Aided Design & Computer Graphics, 2003, 15(1): 1-13.

[4] CUILLIERE J C.A direct method for the automatic discretization of 3D parametric curves [J].Computer Aided Design, 1997, 29(9): 639-647.

[5] CUILLIERE J C.An adaptive method for the automatic triangulation of 3D parametric surfaces [J].Computer Aided Design, 1998, 30(2): 139-149.

[6] LEE C K, HOBBS R E.Automatic adaptive finite element mesh generation over arbitrary two-dimensional domain using advancing front technique [J].Computers & Structures, 1999, 71(1): 9-34.

[7] LEE C K.Automatic metric advancing front triangulation over curved surfaces [J].Engineering Computations, 2000, 17(1): 48-74.

[8] LO S H.Automatic mesh generation over intersecting surfaces [J].International Journal for Numerical Methods in Engineering, 1995, 38(6): 943-954.

[9] ITO Y, NAKAHASHI K.Surface triangulation for polygonal models based on CAD data [J].International Journal for Numerical Methods in Fluids, 2002, 39(1): 75-96.

[10]董亮, 劉厚林, 談明高, 等.離心泵四面體網格質量衡量準則及優化算法 [J].西安交通大學學報, 2011, 45(11): 100-105.DONG Liang, LIU Houlin, TAN Minggao, et al.Quality measurement criteria and optimization algorithm of tetrahedral mesh for centrifugal pumps [J].Journal of Xi’an Jiaotong University, 2011, 45(11): 100-105.

(編輯 苗凌)

ANewAlgorithmforTriangularMeshGenerationbyAdvancingFrontTechnique

WANG Wei,FAN Hongzhou,XI Guang

(School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China)

A new algorithm for triangular mesh generation by advancing front technique is proposed to remove the mapping distortion when the parametric mesh is generated by mapping method.For the current node, this algorithm combines direct method with mapping method to directly locate a new node and to determine the validity in the real space, the mesh in the parametric spaces is generated by mapping method.The angle formed by the adjacent fronts is divided into several units according to the angle scope.The acute angles of nodes is firstly divided, to always retain the polygon fronts with obtuse angle.The proposed algorithm enables to reduce the content and number of validities, avoid the repeated subdivision, cancels the part of the convergence checking, and generate a highly quality mesh.The mesh arrangement can also be used to determine spatial curvature distribution of 3D surfaces for analyzing vibration and precision cutting of the blade.

finite element method; surface mesh; advancing front technique; triangulation; mapping distortion

10.7652/xjtuxb201403012

2013-07-09。

王偉(1982—),男,博士生;席光(通信作者),男,教授,博士生導師。

國家“973計劃”資助項目(2011CB706505);陜西省自然科學基金資助項目(2012JM7004)。

時間: 2013-12-19

TP391

:A

:0253-987X(2014)03-0061-07

網絡出版地址: http:∥www.cnki.net/kcms/detail/61.1069.T.20131219.1121.003.html