徐 青，奚鵬飛

(武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072)

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

網(wǎng)格生成方法有結(jié)構(gòu)化網(wǎng)格生成方法和非結(jié)構(gòu)化網(wǎng)格生成方法。結(jié)構(gòu)化網(wǎng)格生成方法主要有坐標轉(zhuǎn)換法、Paving法等，生成的網(wǎng)格質(zhì)量好，但僅適用于簡單域，一般用在流體動力學(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]。對于自由曲面網(wǎng)格的劃分與優(yōu)化，也有學(xué)者開展了專門研究[3,4]。

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

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

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

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

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

1.1 區(qū)域描述

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

1.2 背景網(wǎng)格

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

1.3 邊界離散

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

(1)

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

1.4 節(jié)理處理

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

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

1.5 網(wǎng)格生成

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

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

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

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

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

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

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

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

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

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

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

(2)

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

1.7.2 拓撲關(guān)系優(yōu)化

當網(wǎng)格局部區(qū)域有一些節(jié)點的度值不理想時，無法僅通過移動節(jié)點的位置來提髙網(wǎng)格的質(zhì)量，需要先改善網(wǎng)格的拓撲結(jié)構(gòu)，才能更好地提髙網(wǎng)格的質(zhì)量。

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

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

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

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

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

圖3 拓撲優(yōu)化示意圖Fig.3 Pattern of topological optimization

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

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

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

2.1.1 網(wǎng)格嵌入

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

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

圖4 新單元點嵌入相鄰單元中Fig.4 New node inserted into adjacent element

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

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

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

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

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

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

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

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

(1)對節(jié)理邊不進行延伸。將節(jié)理邊細分，并生成節(jié)理單元，同時更新節(jié)理邊相鄰的邊端節(jié)點。

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

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

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

圖7 補充邊示意圖Fig.7 Supplementary side

2.2 網(wǎng)格質(zhì)量優(yōu)化改進

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

如圖8所示，由圖8(a)到圖8(b)，表示拓撲關(guān)系優(yōu)化中的節(jié)點刪除，但由于該操作又產(chǎn)生了新的低質(zhì)量單元，且單元中存在邊界節(jié)點，拓撲關(guān)系優(yōu)化操作受到限制。改進思路是：從單元的內(nèi)部節(jié)點入手，調(diào)整節(jié)點A至節(jié)點C，將該單元直接刪除，如圖8(c)所示，周圍單元拓撲關(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個面域，80條邊，65個端點。各區(qū)域材料的滲透系數(shù)如表1所示。由于逸出點處及混凝土防滲墻處滲透坡降較大，表1中也給出了相應(yīng)材料的允許滲透坡降，用以判斷土石壩的滲透穩(wěn)定性。

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

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

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

表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é) 語

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

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

□