四邊形網(wǎng)格自動(dòng)生成方法改進(jìn)及工程應(yīng)用

徐 青，奚鵬飛

(武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072)

隨著計(jì)算機(jī)技術(shù)和數(shù)值分析理論的發(fā)展，有限單元法在工程結(jié)構(gòu)分析中的應(yīng)用越來(lái)越廣泛。有限元網(wǎng)格劃分是有限元分析前處理的主要工作內(nèi)容，其復(fù)雜性及魯棒性一直以來(lái)都是有限元分析的瓶頸問(wèn)題，一定程度上影響著有限單元法的推廣應(yīng)用和發(fā)展[1]。

網(wǎng)格生成方法有結(jié)構(gòu)化網(wǎng)格生成方法和非結(jié)構(gòu)化網(wǎng)格生成方法。結(jié)構(gòu)化網(wǎng)格生成方法主要有坐標(biāo)轉(zhuǎn)換法、Paving法等，生成的網(wǎng)格質(zhì)量好，但僅適用于簡(jiǎn)單域，一般用在流體動(dòng)力學(xué)分析中；非結(jié)構(gòu)化網(wǎng)格生成方法主要有八叉樹法(Octree method)、Delauney法、Octree/Delauney法以及行波法(advancing front method)等，生成的網(wǎng)格質(zhì)量要求寬松，可適應(yīng)復(fù)雜域，主要用于結(jié)構(gòu)數(shù)值分析中[2]。對(duì)于自由曲面網(wǎng)格的劃分與優(yōu)化，也有學(xué)者開展了專門研究[3,4]。

坐標(biāo)轉(zhuǎn)換法[5]是通過(guò)坐標(biāo)轉(zhuǎn)換，將不規(guī)則域轉(zhuǎn)換成規(guī)則域，在規(guī)則域上進(jìn)行網(wǎng)格生成，然后轉(zhuǎn)換回不規(guī)則區(qū)域。Paving法[6]是由邊界向內(nèi)部逐排生成單元，對(duì)邊界的適應(yīng)性較強(qiáng)，但域內(nèi)單元密度由邊界決定，難以生成梯度變化較大的網(wǎng)格。八叉樹法[7]是先將研究域劃分為子域，實(shí)現(xiàn)對(duì)復(fù)雜邊界的分解，再按照一定的要求劃分成單元。該法可以實(shí)現(xiàn)對(duì)網(wǎng)格生成進(jìn)程的控制，邊界適應(yīng)性強(qiáng)，離散區(qū)域的幾何拓?fù)潢P(guān)系不變，但八叉樹法僅僅是一種邊界簡(jiǎn)化方法，不是一種真正的網(wǎng)格生成方法，須依賴于其他的網(wǎng)格生成方法。Delauney法[8]概念簡(jiǎn)單，易于實(shí)施，可以得到網(wǎng)格尺度場(chǎng)梯度變化較大的網(wǎng)格，但邊界和域內(nèi)會(huì)生成可能引起計(jì)算問(wèn)題的點(diǎn)，以及得不到所需網(wǎng)格密度等問(wèn)題。1987年，Peraire等人提出了行波法[9]。行波法是在生成結(jié)點(diǎn)的同時(shí)生成單元網(wǎng)格，應(yīng)用背景網(wǎng)格存放離散所需要的信息，能夠很好地控制所生成的四面體單元的形狀。

陳勝宏等在進(jìn)行系統(tǒng)的自適應(yīng)有限單元分析時(shí)，采用行波法實(shí)現(xiàn)了二維四邊形、三維四面體及三維六面體網(wǎng)格自動(dòng)離散[2,10-12]，編制了core-AFEM系列軟件。本文基于已有研究成果，針對(duì)水工建筑物中具有水平邊和鉛直邊的薄層構(gòu)造及其他復(fù)雜邊界，在進(jìn)行網(wǎng)格自動(dòng)劃分時(shí)可能出現(xiàn)異常中斷等問(wèn)題，探討二維四邊形網(wǎng)格自動(dòng)生成的優(yōu)化及改進(jìn)方法，提高程序的魯棒性，優(yōu)化網(wǎng)格質(zhì)量。

1 平面四邊形網(wǎng)格生成

水工結(jié)構(gòu)及巖土工程，其研究區(qū)域通常由多種材料構(gòu)成，且存在斷層、軟弱夾層、防滲墻等薄層結(jié)構(gòu)，以及廊道(孔洞)等，使得網(wǎng)格劃分形狀不規(guī)則，隨機(jī)性較大。針對(duì)這些特點(diǎn)，網(wǎng)格生成的總體思想是，將多介質(zhì)復(fù)雜區(qū)域轉(zhuǎn)換為若干由單一介質(zhì)構(gòu)成的子域，每一子域單獨(dú)生成網(wǎng)格，由各子域的網(wǎng)格組合成復(fù)雜區(qū)域的網(wǎng)格。

圖1為四邊形網(wǎng)格劃分流程，主要技術(shù)特點(diǎn)是包含節(jié)理邊處理、邊界離散、三角單元生成過(guò)程中節(jié)點(diǎn)選取、三角形合并為四邊形的控制過(guò)程，以及網(wǎng)格質(zhì)量的初步優(yōu)化。

1.1 區(qū)域描述

計(jì)算域由子域、邊界和節(jié)點(diǎn)3級(jí)數(shù)據(jù)結(jié)構(gòu)進(jìn)行描述。首先根據(jù)材料分區(qū)和結(jié)構(gòu)分區(qū)，將復(fù)雜區(qū)域分解成若干個(gè)簡(jiǎn)單的子區(qū)域，由邊界線描述子區(qū)域，再由節(jié)點(diǎn)構(gòu)造邊界線，最后把由邊界線描述的子區(qū)域轉(zhuǎn)換為由其邊界節(jié)點(diǎn)描述。邊界節(jié)點(diǎn)按逆時(shí)針?lè)较蚺帕校瑓^(qū)域位于節(jié)點(diǎn)排列方向的左側(cè)，逆向排列的節(jié)點(diǎn)封閉環(huán)唯一確定了區(qū)域。

1.2 背景網(wǎng)格

背景網(wǎng)格是一套覆蓋全部區(qū)域的尺度網(wǎng)格，用于給定區(qū)域內(nèi)單元尺寸的分布信息。網(wǎng)格上每一節(jié)點(diǎn)給定了該點(diǎn)的相鄰單元尺寸，區(qū)域內(nèi)任一點(diǎn)的相鄰單元尺寸由該點(diǎn)所在的背景網(wǎng)格單元插值確定。

1.3 邊界離散

邊界線準(zhǔn)備數(shù)據(jù)文件中，邊界點(diǎn)是區(qū)域幾何形態(tài)等發(fā)生變化的點(diǎn)，一般不直接用于網(wǎng)格生成，而是依據(jù)背景網(wǎng)格，先離散邊界線，增加節(jié)點(diǎn)數(shù)。對(duì)于任意線段，其上任意一點(diǎn)的尺度可由該點(diǎn)所在的背景單元內(nèi)插確定。設(shè)δ(l)為單元尺寸沿邊界線l的分布函數(shù)，L為l的長(zhǎng)度，邊界節(jié)點(diǎn)加密后邊界線l的線段數(shù)N取與N′最接近的整數(shù)，N′由下式計(jì)算得到：

(1)

子域內(nèi)生成全四邊形網(wǎng)格的必要條件是：邊界節(jié)點(diǎn)數(shù)必須是偶數(shù)。對(duì)于任一子區(qū)域，選取一個(gè)最長(zhǎng)待劃分邊，通常是劃分份數(shù)最多的普通邊，該劃分邊的劃分份數(shù)可能會(huì)增減1，以保證總節(jié)點(diǎn)數(shù)為偶數(shù)。

1.4 節(jié)理處理

研究域中存在節(jié)理、防滲墻等薄層結(jié)構(gòu)時(shí)，統(tǒng)一定義其中線為節(jié)理邊。先對(duì)邊界進(jìn)行離散，再根據(jù)薄層結(jié)構(gòu)的厚度t，生成節(jié)理單元，如圖2所示。

圖2 節(jié)理單元生成Fig.2 Joint element generation

1.5 網(wǎng)格生成

網(wǎng)格自動(dòng)劃分算法以三角形單元生成為基礎(chǔ)，設(shè)置四邊形單元的狀態(tài)參數(shù)，在一個(gè)三角形單元基礎(chǔ)上生成相鄰三角形單元，從而構(gòu)成四邊形單元，完成全域的四邊形網(wǎng)格劃分。

新生成的單元需進(jìn)行單元質(zhì)量控制，主要包括有效性檢驗(yàn)、單元形狀控制和奇異性檢驗(yàn)。若有一項(xiàng)為否，則跳到下一節(jié)點(diǎn)再次判斷。若全部滿足，則該節(jié)點(diǎn)與活動(dòng)邊形成新的三角形單元。

1.6 網(wǎng)格質(zhì)量控制

(1)有效性控制。網(wǎng)格生成過(guò)程中，需對(duì)其有效性進(jìn)行控制，新生成的三角形網(wǎng)格不能穿過(guò)已生成的四邊形網(wǎng)格，以保證網(wǎng)格的合理性和有效性。

(2)四邊形內(nèi)角控制。四邊形單元的最優(yōu)形式是正方形，因此控制四邊形的內(nèi)角可以對(duì)網(wǎng)格的質(zhì)量進(jìn)行基本控制，本程序中控制四邊形的各內(nèi)角不大于120°。若不滿足，則判定網(wǎng)格質(zhì)量為不合格。

(3)奇異單元控制。對(duì)于邊界單元，尤其是在邊相鄰的區(qū)域單元，可能會(huì)出現(xiàn)三邊四結(jié)點(diǎn)的三角形奇異單元，對(duì)此本文通過(guò)控制四邊形單元中各邊節(jié)點(diǎn)數(shù)不超過(guò)2來(lái)避免此類單元的產(chǎn)生。

1.7 網(wǎng)格質(zhì)量?jī)?yōu)化

對(duì)于復(fù)雜域，通常會(huì)出現(xiàn)自動(dòng)生成的網(wǎng)格質(zhì)量不高的現(xiàn)象，因此需對(duì)網(wǎng)格進(jìn)行優(yōu)化。優(yōu)化方式包括節(jié)點(diǎn)位置優(yōu)化和拓?fù)潢P(guān)系優(yōu)化2種。

1.7.1 節(jié)點(diǎn)位置優(yōu)化

節(jié)點(diǎn)位置優(yōu)化是通過(guò)保持網(wǎng)格拓?fù)潢P(guān)系不變的情況下改變節(jié)點(diǎn)的位置，以此來(lái)提高網(wǎng)格質(zhì)量的一種技術(shù)。其中最經(jīng)典的是拉普拉斯光滑算法，該算法簡(jiǎn)單，且耗時(shí)短。

節(jié)點(diǎn)位置優(yōu)化通過(guò)調(diào)整網(wǎng)格單元內(nèi)部的節(jié)點(diǎn)或者表面非邊界節(jié)點(diǎn)實(shí)現(xiàn)。假設(shè)要移動(dòng)的節(jié)點(diǎn)為K，優(yōu)化后的K點(diǎn)坐標(biāo)為：

(2)

式中：Ni為該節(jié)點(diǎn)所在單元內(nèi)的其他節(jié)點(diǎn)總數(shù)。

1.7.2 拓?fù)潢P(guān)系優(yōu)化

當(dāng)網(wǎng)格局部區(qū)域有一些節(jié)點(diǎn)的度值不理想時(shí)，無(wú)法僅通過(guò)移動(dòng)節(jié)點(diǎn)的位置來(lái)提髙網(wǎng)格的質(zhì)量，需要先改善網(wǎng)格的拓?fù)浣Y(jié)構(gòu)，才能更好地提髙網(wǎng)格的質(zhì)量。

網(wǎng)格的邊界節(jié)點(diǎn)是不變的，因此拓?fù)潢P(guān)系優(yōu)化主要針對(duì)內(nèi)部節(jié)點(diǎn)進(jìn)行。對(duì)四邊形網(wǎng)格而言，內(nèi)部節(jié)點(diǎn)的相鄰單元個(gè)數(shù)NE最好等于4。當(dāng)某一個(gè)內(nèi)部節(jié)點(diǎn)的相鄰單元個(gè)數(shù)NE比4大或小很多時(shí)，環(huán)繞該節(jié)點(diǎn)的單元?jiǎng)t會(huì)產(chǎn)生較大的畸變，此時(shí)應(yīng)進(jìn)行調(diào)整。調(diào)整方法主要有(見圖3)。

(1)節(jié)點(diǎn)刪除[圖3(a)]。如果內(nèi)部節(jié)點(diǎn)A的相鄰單元數(shù)NE=2，則節(jié)點(diǎn)A將被刪除。

(2)單元?jiǎng)h除[圖3(b)]。如果某單元的2個(gè)相對(duì)節(jié)點(diǎn)A、B都是內(nèi)部節(jié)點(diǎn)，且NE(A)=NE(B)=3，則刪除該單元。

(3)對(duì)角線交換[圖3(c)]。由于四邊形網(wǎng)格中每個(gè)節(jié)點(diǎn)所在邊的個(gè)數(shù)NE等于4時(shí)，網(wǎng)格質(zhì)量最優(yōu)。因此，首先檢查由2個(gè)內(nèi)部節(jié)點(diǎn)連成的邊AB，若NE(A)+NE(B)≥9，再檢查以AB為邊的2個(gè)單元組成的六邊形AEDBFC中其他對(duì)角線的2個(gè)節(jié)點(diǎn)是否更接近最優(yōu)。有3種情況：當(dāng)[NE(A)+NE(B)]-[NE(C)+NE(D)]>[NE(A)+NE(B)]-[NE(E)+NE(F)]≥3時(shí)，則用CD代替AB；當(dāng)[NE(A)+NE(B)]-[NE(E)+NE(F)]>[NE(A)+NE(B)]-[NE(C)+NE(D)]≥3時(shí)，則用EF代替AB；否則不進(jìn)行對(duì)角線交換。

(4)邊刪除[圖3(d)]。檢查所有2個(gè)內(nèi)部節(jié)點(diǎn)連成的邊EF，若NE(A)=NE(B)=3，且含有EF邊的2個(gè)相鄰單元中沒有邊界邊，則刪除邊EF。當(dāng)NE(A)+NE(D)≤NE(B)+NE(C)時(shí)，用AD作為新的替代邊；否則用BC作為新的替代邊。

圖3 拓?fù)鋬?yōu)化示意圖Fig.3 Pattern of topological optimization

2 四邊形網(wǎng)格自動(dòng)生成算法改進(jìn)及優(yōu)化

在壩工工程中，防滲墻是常用的工程措施，與大壩相比，防滲墻屬于薄層結(jié)構(gòu)；同時(shí)，壩基巖體通常存在著各類軟弱夾層。本文針對(duì)具有水平邊或垂直邊的薄層單元等，提出四邊形網(wǎng)格生成改進(jìn)方法，以提高算法的魯棒性，并進(jìn)一步優(yōu)化網(wǎng)格質(zhì)量。

2.1 四邊形網(wǎng)格生成算法改進(jìn)

2.1.1 網(wǎng)格嵌入

網(wǎng)格嵌入或交叉是網(wǎng)格生成過(guò)程中最為常見的情形，通常是由于算法邏輯沒有考慮到某些特殊情形(如具有水平邊或鉛直邊的薄層單元)，交叉判定時(shí)產(chǎn)生錯(cuò)誤，從而影響網(wǎng)格魯棒生成。

(1)水平薄層單元。如圖4所示，N1、N2為當(dāng)前活動(dòng)波前，(X,Y)為可能的新單元點(diǎn)。首先考察點(diǎn)P是否滿足新點(diǎn)要求。從圖4可以清晰地看到，點(diǎn)P與N1、N2構(gòu)成的新三角嵌入到已生成單元中(圖4中陰影部分)，不滿足網(wǎng)格的有效性要求，故應(yīng)排除該點(diǎn)，并進(jìn)入下一階段新單元點(diǎn)選取。

圖4 新單元點(diǎn)嵌入相鄰單元中Fig.4 New node inserted into adjacent element

引起該嵌入發(fā)生的原因是：線段PM與線段BN1的交叉判定出錯(cuò)。在計(jì)算交點(diǎn)J的過(guò)程中，由于舍入誤差的存在，計(jì)算所得的J點(diǎn)縱坐標(biāo)Y可能出現(xiàn)無(wú)限接近卻不等于點(diǎn)P和點(diǎn)M的縱坐標(biāo)的情況(PM是一水平線)，導(dǎo)致J點(diǎn)被認(rèn)為不在PM上，從而該三角單元被允許生成，形成嵌入。對(duì)于由于計(jì)算誤差產(chǎn)生的線段相交誤判，本文提出的改進(jìn)方法是：對(duì)線段相交模塊進(jìn)行優(yōu)化，在直線相交子程序中，當(dāng)存在水平邊時(shí)，直接用水平邊的縱坐標(biāo)指定交點(diǎn)Y值，再計(jì)算交點(diǎn)X值，從而避免線段相交誤判。

(2)鉛直薄層單元。如圖5所示，N1、N2為當(dāng)前活動(dòng)波前，虛線為節(jié)理邊。點(diǎn)P為可能的新的單元點(diǎn)。考察點(diǎn)P是否滿足新點(diǎn)要求，從圖5可以清晰地看到：P與N1、N2構(gòu)成的新三角嵌入到節(jié)理單元中，不滿足網(wǎng)格的有效性要求。

圖5 新單元點(diǎn)嵌入節(jié)理單元中Fig.5 New node inserted into joint element

引起該嵌入發(fā)生的原因是：線段PM與線段N3N4的交叉判定出錯(cuò)。這里存在2次舍入誤差。第1次舍入誤差是N4的坐標(biāo)計(jì)算，N4是節(jié)理擴(kuò)展邊上的節(jié)點(diǎn)。在計(jì)算該點(diǎn)橫坐標(biāo)的過(guò)程中，由于舍入誤差的存在，使得計(jì)算得到的節(jié)點(diǎn)N4的橫坐標(biāo)也是與N3橫坐標(biāo)十分接近的一個(gè)數(shù)值，因此N3N4線段的節(jié)點(diǎn)斜率絕對(duì)值極大。第2次舍入誤差來(lái)源于PM與線段N3N4計(jì)算所得交點(diǎn)橫坐標(biāo)，十分接近，但是不相等，計(jì)算縱坐標(biāo)時(shí)由于斜率極大，使得縱坐標(biāo)與理論結(jié)果相差較大，得出相交節(jié)點(diǎn)不在N3N4上的結(jié)論，從而生成嵌入單元。

2.1.2 內(nèi)部斷層單元處理優(yōu)化

在有限元數(shù)值模型建模階段，常常將內(nèi)部節(jié)理邊延伸至模型邊界，這樣不可避免地會(huì)影響周圍區(qū)域的波前，使網(wǎng)格劃分更加密集，增加網(wǎng)格數(shù)量，見圖6。

圖6 斷層邊延伸至邊界Fig.6 Extend the joint edge to the boundary

本文對(duì)內(nèi)部節(jié)理邊的處理提出如下優(yōu)化方案。

(1)對(duì)節(jié)理邊不進(jìn)行延伸。將節(jié)理邊細(xì)分，并生成節(jié)理單元，同時(shí)更新節(jié)理邊相鄰的邊端節(jié)點(diǎn)。

(2)對(duì)邊界與節(jié)理邊垂直的面域B和E，檢查面域的連續(xù)性。這種面域在節(jié)理單元生成、節(jié)理邊端點(diǎn)更新后，面域端點(diǎn)更新為節(jié)理單元兩端，整個(gè)面域波前不再連續(xù)。

(3)增加補(bǔ)充邊。增加邊12和邊34，如圖7所示。將該補(bǔ)充邊默認(rèn)為已劃分狀態(tài)，避免被后續(xù)可能的邊界加密，產(chǎn)生額外的節(jié)點(diǎn)，破壞已生成的節(jié)理單元。

(4)為避免補(bǔ)充邊的生成引起受影響面域的波前總節(jié)點(diǎn)數(shù)由偶數(shù)變?yōu)槠鏀?shù)，四邊形網(wǎng)格自動(dòng)生成的前提條件不滿足，使得網(wǎng)格劃分失敗，須將操作過(guò)程與邊界離散操作的順序進(jìn)行交換。

圖7 補(bǔ)充邊示意圖Fig.7 Supplementary side

2.2 網(wǎng)格質(zhì)量?jī)?yōu)化改進(jìn)

在網(wǎng)格質(zhì)量?jī)?yōu)化的過(guò)程中，網(wǎng)格拓?fù)潢P(guān)系發(fā)生改變，從而影響周圍網(wǎng)格的拓?fù)潢P(guān)系，可能會(huì)產(chǎn)生新的不良單元。對(duì)此，本文提出了相應(yīng)的優(yōu)化方案。

如圖8所示，由圖8(a)到圖8(b)，表示拓?fù)潢P(guān)系優(yōu)化中的節(jié)點(diǎn)刪除，但由于該操作又產(chǎn)生了新的低質(zhì)量單元，且單元中存在邊界節(jié)點(diǎn)，拓?fù)潢P(guān)系優(yōu)化操作受到限制。改進(jìn)思路是：從單元的內(nèi)部節(jié)點(diǎn)入手，調(diào)整節(jié)點(diǎn)A至節(jié)點(diǎn)C，將該單元直接刪除，如圖8(c)所示，周圍單元拓?fù)潢P(guān)系不變，網(wǎng)格質(zhì)量得到優(yōu)化。

圖8 邊界三角單元處理Fig.8 Treatment of singular element on boundary

3 工程應(yīng)用

本文選取具有復(fù)雜地質(zhì)構(gòu)造，并設(shè)置有混凝土防滲墻的土石壩剖面，如圖9所示。將其分為16個(gè)面域，80條邊，65個(gè)端點(diǎn)。各區(qū)域材料的滲透系數(shù)如表1所示。由于逸出點(diǎn)處及混凝土防滲墻處滲透坡降較大，表1中也給出了相應(yīng)材料的允許滲透坡降，用以判斷土石壩的滲透穩(wěn)定性。

圖9 土石壩剖面(單位：m)Fig.9 Section of an embankment dam

依據(jù)上述改進(jìn)方法，四邊形網(wǎng)格自動(dòng)劃分過(guò)程順利，網(wǎng)格質(zhì)量良好，如圖10所示。

采用課題組開發(fā)的軟件core-AFEM/seep2對(duì)該土石壩進(jìn)行滲流分析，壩體自由面、等水頭線及滲透坡降、流速如圖11～圖13所示。計(jì)算結(jié)果顯示：①自由面及各滲透特性分布規(guī)律符合常規(guī)；②J混凝土防滲墻≤40，J二期填土<0.48，vmax<10-6m/s，說(shuō)明滿足滲透穩(wěn)定要求，與實(shí)際工程狀況相符。

表1 各區(qū)域材料滲透系數(shù)Tab.1 The permeability coefficient

圖10 四邊形網(wǎng)格劃分結(jié)果Fig.10 The quadrilateral mesh

圖11 等水頭線分布Fig.11 The equipotential line distribution

圖12 滲透坡降分布Fig.12 The seepage gradient distribution

4 結(jié) 語(yǔ)

本文基于行波法，針對(duì)具有水平邊和鉛直邊的薄層單元(如混凝土防滲墻、斷層等)，提出了二維四邊形網(wǎng)格自動(dòng)劃分改進(jìn)方法，提高了網(wǎng)格自動(dòng)生成的魯棒性，并經(jīng)工程應(yīng)用檢驗(yàn)，可以很好地應(yīng)用于有限元數(shù)值分析，并具有網(wǎng)格劃分合理、易于掌握等特點(diǎn)，極大地提高了網(wǎng)格自動(dòng)生成的適用范圍。

圖13 流速分布Fig.13 The flow velocity distribution

□