多介質(zhì)復(fù)雜區(qū)域四邊形網(wǎng)格自動(dòng)剖分算法及應(yīng)用

付成華，周洪波

(1.西華大學(xué)能源與環(huán)境學(xué)院，成都 610039;2.二灘水電開(kāi)發(fā)有限責(zé)任公司，成都 610051)

1 概述

有限元計(jì)算程序不斷涌現(xiàn)出來(lái)并成功應(yīng)用于實(shí)際工程中，圍繞有限元計(jì)算程序集成開(kāi)發(fā)的軟件也越來(lái)越多。在有限元分析中，網(wǎng)格生成技術(shù)一直是制約科學(xué)和工程計(jì)算精度和效率的一個(gè)重要因素，所以根據(jù)計(jì)算區(qū)域的幾何結(jié)構(gòu)信息，利用計(jì)算機(jī)進(jìn)行自動(dòng)剖分，產(chǎn)生有限元計(jì)算所必須的網(wǎng)格數(shù)據(jù)，進(jìn)一步提高有限元數(shù)值計(jì)算的效率，加大有限單元法及其集成軟件的推廣應(yīng)用是非常必要的。

本文提出一種四邊形網(wǎng)格自動(dòng)剖分的改進(jìn)行波算法。無(wú)需設(shè)置背景網(wǎng)格，將多介質(zhì)復(fù)雜區(qū)域分成多個(gè)封閉子域，先生成平面三角形網(wǎng)格，優(yōu)化后通過(guò)合并、分解生成四邊形單元，再次對(duì)四邊形單元優(yōu)化處理即可獲得過(guò)渡光滑且質(zhì)量良好的四邊形網(wǎng)格，從而實(shí)現(xiàn)多介質(zhì)復(fù)雜區(qū)域的四邊形網(wǎng)格自動(dòng)剖分。最后通過(guò)應(yīng)用于某鐵路路堤斷面和瀑布溝電站廠房2#機(jī)組斷面計(jì)算域的網(wǎng)格剖分，驗(yàn)證了該方法的可行性及可靠性。

2 算法步驟

行波法最初由Lo提出并用于平面三角形網(wǎng)格的自動(dòng)生成，后來(lái)通過(guò)引入背景網(wǎng)格來(lái)控制內(nèi)部節(jié)點(diǎn)的生成，并推廣到任意形狀三維區(qū)域的四面體網(wǎng)格生成。行波法主要從區(qū)域邊界開(kāi)始，定義一個(gè)“前沿”，在前沿上滿足一定條件的地方生成一系列網(wǎng)格，同時(shí)更新“前沿”，如此不斷地重復(fù)循環(huán)，直到前沿縮為一點(diǎn)，結(jié)束網(wǎng)格的生成。行波法剖分網(wǎng)格的自動(dòng)化程度高，對(duì)復(fù)雜幾何形狀與邊界的二維(或三維)區(qū)域中的三角形(或四面體)網(wǎng)格生成具有很大的靈活性及可靠性［1－2］。

本文改進(jìn)行波法是基于計(jì)算域AutoCAD模型和程序語(yǔ)言實(shí)現(xiàn)的，四邊形網(wǎng)格自動(dòng)剖分的具體步驟如下:

(1)將多介質(zhì)復(fù)雜區(qū)域看成是多個(gè)封閉子區(qū)域的組合，利用編制的AutoCAD應(yīng)用程序ptline.lsp從計(jì)算域的AutoCAD圖中提取各封閉子區(qū)域邊界的點(diǎn)線信息。

(2)離散區(qū)域邊界，以離散的各子域邊界作為“前沿”，利用改進(jìn)行波法在各個(gè)子域內(nèi)生成三角形網(wǎng)格，并進(jìn)行三角形網(wǎng)格質(zhì)量評(píng)價(jià)和優(yōu)化。

(3)分子域?qū)?個(gè)三角形合并為1個(gè)四邊形，合并后的四邊形分解為2個(gè)或4個(gè)小四邊形，未合并三角形分解為3個(gè)四邊形。

(4)對(duì)四邊形網(wǎng)格進(jìn)行質(zhì)量評(píng)價(jià)及優(yōu)化，實(shí)現(xiàn)整個(gè)計(jì)算域四邊形網(wǎng)格的自動(dòng)剖分。

3 關(guān)鍵技術(shù)改進(jìn)

3.1 子域內(nèi)部節(jié)點(diǎn)生成

將多介質(zhì)復(fù)雜區(qū)域看作多個(gè)封閉子域進(jìn)行剖分，每個(gè)子域在剖分過(guò)程中又被分為已剖分域和未剖分域，封閉域內(nèi)部節(jié)點(diǎn)的生成步驟(圖1)如下:

(1)根據(jù)網(wǎng)格密度控制參數(shù)離散各子域邊界線得到各子域的初始生成波。

(2)從子域生成波中選取內(nèi)角最小的點(diǎn)，作為新三角形的基點(diǎn)。

若基點(diǎn)A內(nèi)角小于等于90°，則與基點(diǎn)相鄰的2點(diǎn)B，C直接作為新三角形的頂點(diǎn)，生成三角形ABC，記錄子域未剖分域，更新子域生成波。

若基點(diǎn)B內(nèi)角＞90°，在∠DBC平分線上取點(diǎn)P，LBP=(LBD+LBC)/2，且P點(diǎn)在子域內(nèi)部。若P點(diǎn)不在子域內(nèi)，則取∠DBC平分線與子域邊界的交點(diǎn)為P。判別∠FDP和∠PCE，若∠FDP、∠PCE均≥30°，則P點(diǎn)作為新生成的內(nèi)部節(jié)點(diǎn)，生成三角形 BDP 和BPC;若∠FDP＜30°、∠PCE≥30°，則直接生成三角形BDF和BFC;若∠FDP≥30°、∠PCE＜30°，則直接生成三角形BDE和BEC。

(3)記錄子域未剖分域，更新子域生成波。重復(fù)(2)，直至生成波為空。

圖1 內(nèi)部節(jié)點(diǎn)生成示意圖Fig.1 Schematic diagram of the generation of internal nodes

3.2 網(wǎng)格密度控制

網(wǎng)格密度指單位面積內(nèi)單元的平均個(gè)數(shù)，這里網(wǎng)格密度控制采用總體密度控制和局部過(guò)度控制相結(jié)合的方式實(shí)現(xiàn)。

(1)總體密度控制。由用戶根據(jù)自身需要設(shè)定待剖分區(qū)域網(wǎng)格總體密度的控制參數(shù)，參數(shù)在(0，1)之間取值，數(shù)值越小說(shuō)明網(wǎng)格越密。

(2)局部過(guò)渡控制。各子域按面積由小到大排序，依次離散各子域的邊界線。如圖2，離散子域邊界線AB時(shí)，以2端點(diǎn)A，B為垂足分別向子域內(nèi)做垂線AI和BJ，與子域邊界的交點(diǎn)為I和J。A端離散長(zhǎng)度l取總體控制長(zhǎng)度、相鄰邊界線FA的長(zhǎng)度 LFA和LAI三者中的最小者。

圖2 邊界離散控制示意圖Fig.2 Schematic diagram of border discrete control

3.3 網(wǎng)格質(zhì)量評(píng)價(jià)及優(yōu)化

網(wǎng)格自動(dòng)剖分過(guò)程中，單元質(zhì)量控制非常重要。這里通過(guò)建立網(wǎng)格質(zhì)量評(píng)價(jià)體系和多重優(yōu)化措施來(lái)保證最終獲得良好的四邊形網(wǎng)格。

3.3.1 網(wǎng)格質(zhì)量評(píng)價(jià)

(1)三角形網(wǎng)格質(zhì)量評(píng)價(jià)。三角形單元質(zhì)量主要通過(guò)單元內(nèi)角來(lái)判斷。三角形單元內(nèi)角介于45°和135°之間，質(zhì)量好;內(nèi)角都介于30°和150°之間，質(zhì)量一般;若最小內(nèi)角＜30°或最大內(nèi)角＞150°，質(zhì)量較差［1］。

(2)四邊形網(wǎng)格質(zhì)量評(píng)價(jià)。按逆時(shí)針走向的四邊形ABCD被對(duì)角線AC分成2個(gè)三角形ABC和ACD，由公式(1)計(jì)算其 α值分別為和;同時(shí)四邊形ABCD又可被對(duì)角線BD分成2個(gè)三角形ABD和BCD，其α值分別為和。將4個(gè)三角形的α 值，，按降序重新排列成為α1，α2，α3，α4。則四邊形 ABCD的質(zhì)量判斷采用文獻(xiàn)［3－4］準(zhǔn)則:

矩形的β值為1，β值越小，四邊形質(zhì)量越差。

3.3.2 網(wǎng)格優(yōu)化

(1)三角形網(wǎng)格優(yōu)化。本算法中三角形單元優(yōu)化采用文獻(xiàn)［1］中的優(yōu)化方式:①刪除節(jié)點(diǎn);②對(duì)角線調(diào)整;③拉普拉斯優(yōu)化。

(2)四邊形網(wǎng)格優(yōu)化。對(duì)網(wǎng)格質(zhì)量較差的四邊形采用拉普拉斯光順和內(nèi)部調(diào)整進(jìn)行優(yōu)化［1］。

4 應(yīng)用實(shí)例

(1)某鐵路路堤斷面計(jì)算域如圖3所示。包括＜0＞、＜4＞、＜6＞、＜13－1W4＞、＜13－W4＞、＜13－W2＞ 共6個(gè)不規(guī)則地層，P1，P2和P3，P4之間分別為2軌道寬度。加載應(yīng)用程序ptline.lsp從AutoCAD模型提取6個(gè)封閉子域的點(diǎn)線信息，運(yùn)行自動(dòng)剖分程序即得計(jì)算域的網(wǎng)格自動(dòng)剖分結(jié)果如圖4至圖6所示，粗網(wǎng)格四邊形單元數(shù)215，節(jié)點(diǎn)數(shù)249;細(xì)網(wǎng)格四邊形單元數(shù)756，節(jié)點(diǎn)數(shù)811。優(yōu)化前后單元質(zhì)量分布情況見(jiàn)表1。

(2)瀑布溝電站地下廠房2#機(jī)組剖面計(jì)算域如圖7至圖8所示:水平方向沿主廠房中心線向上游取150m;向下游取250m;鉛直方向由500m高程處取至地表。計(jì)算模型中考慮的斷層、擠壓破碎帶主要包括:f13，f29，f19，g(1)，f(7)，f(8)，f(1)，f28。網(wǎng)格剖分結(jié)果如圖9至圖10所示。優(yōu)化后三角形網(wǎng)格單元數(shù)4 744，節(jié)點(diǎn)數(shù)2 492，其中質(zhì)量好的單元4 382個(gè)，質(zhì)量較差的單元66個(gè)。優(yōu)化后四邊形網(wǎng)格單元數(shù)9 678，節(jié)點(diǎn)數(shù)9 920，其中β≥0.5的單元6 235個(gè)，網(wǎng)格整體剖分效果較好。

圖3 某鐵路路堤斷面計(jì)算域Fig.3 Computational domain of the railway embankment section

圖4 粗網(wǎng)格Fig.4 Coarse grids

圖5 細(xì)網(wǎng)格Fig.5 Fine grids

圖6 優(yōu)化后的局部四角形網(wǎng)格Fig.6 Optimized local quadrilateral grids

表1 優(yōu)化前后單元質(zhì)量分布Table 1 Distribution of quadrilateral element quality before and after the optimization

圖7 瀑布溝電站廠房2#機(jī)組剖面Fig.7 Computatonal domain of the profile of No.2 unit of Pubugou power station underground plant

圖8 四邊形網(wǎng)格Fig.8 Quadrilateral grids

圖9 局部三角形網(wǎng)格Fig.9 Local triangular grids

圖10 局部四邊形網(wǎng)格Fig.10 Local quadrilateral grids

5 結(jié)語(yǔ)

本文提出的四邊形網(wǎng)格自動(dòng)剖分的改進(jìn)行波法，無(wú)需設(shè)置背景網(wǎng)格，直接由AutoCAD模型出發(fā)，加強(qiáng)了網(wǎng)格剖分的過(guò)程控制，能迅速而有力地生成高質(zhì)量的三角形網(wǎng)格和四邊形網(wǎng)格，實(shí)現(xiàn)多介質(zhì)復(fù)雜區(qū)域網(wǎng)格的自動(dòng)剖分，簡(jiǎn)單易行，適應(yīng)性強(qiáng)，既可滿足有限元計(jì)算的需要，又可用于大型有限元計(jì)算系統(tǒng)的集成開(kāi)發(fā)，具有很好的應(yīng)用價(jià)值。

［1］馮道雨.邊坡穩(wěn)定分析的自適應(yīng)有限元法研究［D］.武漢:武漢大學(xué)，2004.(FENG Dao-yu.Self-adaptive Finite Element for Slope Stability Analysis［D］.Wuhan:Wuhan University，2004.(in Chinese))

［2］馮道雨，陳尚法，陳勝宏.復(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.(in Chinese))

［3］劉春太，楊曉東，陳靜波，等，任意平面域自漸變?nèi)切尉W(wǎng)格的自動(dòng)劃分［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，1999，11(4):293－295.(LIU Chun-tai，YANG Xiao-dong，CHEN Jing-bo，etal.Automatic Generation of Self-changing Triangle Grid for Arbitrary Planar Domains［J］.Journal of Computer-Aided Design ＆ Computer Graphics，1999，11(4):293－295.(in Chinese))

［4］楊曉東，申長(zhǎng)雨，陳靜波，等，任意平面區(qū)域的變密度四邊形網(wǎng)格生成方法［J］.工程力學(xué)，2001，18(2):135－139.(YANG Xiao-dong，SHEN Chang-yu，CHEN Jing-bo，etal.Generation of Graded Quadrilateral Meshes for Arbitrary Planar Domains［J］.Engineering Mechanics，2001，18(2):135－139.(in Chinese))

［5］CHEN Jian-jun，ZHENG Yao.A Robust and Quality Guaranteed Pattern Module Scheme for Multi-subdomain Methods in Mesh Generation［J］.Journal of Computer-Aided Design and Computer Graphics，2005，17(10):2286－2292.