基于非線性彈簧近似的平面三角網(wǎng)格優(yōu)化方法

尹 暢，張思全，劉 雨，齊 川

（上海海事大學(xué)物流工程學(xué)院，上海201306）

0 引 言

科學(xué)技術(shù)領(lǐng)域中系統(tǒng)數(shù)學(xué)模型往往采用偏微分方程形式來描述復(fù)雜的物理現(xiàn)象，而數(shù)值計(jì)算技術(shù)是這類數(shù)學(xué)模型的主要計(jì)算分析途徑。PDE數(shù)值解法的核心是將全局定義域中連續(xù)未知函數(shù)轉(zhuǎn)化為一定數(shù)量子定義域上的簡單函數(shù)，進(jìn)一步建立僅僅在定義域內(nèi)一定數(shù)量的離散結(jié)點(diǎn)上處理的代數(shù)方程組。網(wǎng)格劃分（Meshing Generation）就是將PDE全局定義域劃分成離散單元（elements）或者控制體積（control volume）集合，其構(gòu)成了系統(tǒng)的自定義域和離散節(jié)點(diǎn)集合。由于單元采用簡單集合形狀，因此可以方便地在這些子域上構(gòu)造簡單初等函數(shù)逼近連續(xù)的PDE系統(tǒng)響應(yīng)。而且，可以通過調(diào)節(jié)網(wǎng)格類型、單元形狀、尺寸和數(shù)量等因素來提高PDE系統(tǒng)的數(shù)值計(jì)算精度［1］。

三角形是二維空間的單純形，故在二維空間的離散化封面占據(jù)主導(dǎo)地位，相繼出現(xiàn)了插點(diǎn)增量式Delaunay三角形劃分方法、推進(jìn)波前法、分治算法以及掃描算法等。

由于Delaunay三角劃分算法可以使全部三角形的最小內(nèi)角為最大，且保證三角單元不會(huì)重疊（即所有邊最多只被兩個(gè)三角形共有），同時(shí)三角單元的三邊之上不會(huì)有其他節(jié)點(diǎn)，故三角剖分以Delaunay法的應(yīng)用最為廣泛［2］。

1 算法思想與實(shí)現(xiàn)步驟

對形狀規(guī)則的幾何域進(jìn)行網(wǎng)格劃分比較容易，且易于控制單元網(wǎng)格形狀與尺寸，故先建立一個(gè)足夠大的矩形來包圍不規(guī)則幾何域。對矩形使用正交柵格得到節(jié)點(diǎn)p（x，y），利用函數(shù)d（x，y）表示點(diǎn)p與幾何域的位置關(guān)系。在幾何域內(nèi)部及邊界附近的節(jié)點(diǎn)予以保留，與幾何域的頂點(diǎn)和限定點(diǎn)一起作為Delaunay三角劃分的初始點(diǎn)集。

對于不夠飽滿的單元形狀通過彈簧比擬對節(jié)點(diǎn)進(jìn)行受力分析的迭代計(jì)算來使節(jié)點(diǎn)達(dá)到受力平衡，從而提高單元網(wǎng)格的質(zhì)量。將網(wǎng)格節(jié)點(diǎn)視為受力點(diǎn)Pi（i＝1，2，3，．．．，n），每個(gè)Pi都有對應(yīng)的力－位移函數(shù)。受力不平衡（連接到節(jié)點(diǎn)的邊的長度不等）使節(jié)點(diǎn)產(chǎn)生相應(yīng)的位移，同時(shí)利用MATLAB中的Delaunay三角形函數(shù)調(diào)整拓?fù)浣Y(jié)構(gòu)。最終邊長應(yīng)接近所要求的相對長度（相對長度為1時(shí)，則各邊長度幾乎相等）。

1．1 初始離散點(diǎn)集生成

首先建立幾何域的包絡(luò)矩形。然后根據(jù)網(wǎng)格控制尺寸h在縱橫方向布置柵格線，設(shè)置橫向x軸網(wǎng)格尺寸為h，縱向y軸網(wǎng)格尺寸為，然后將偶數(shù)行的點(diǎn)向x軸正向平移h／2，使得相鄰三點(diǎn)構(gòu)成的三角形為理想的正三角形。然后根據(jù)射線法來判定點(diǎn)與幾何域的內(nèi)外關(guān)系［3］，并刪除幾何域外部的全部柵格點(diǎn)，以此來建立Delaunay三角形網(wǎng)格劃分所需的初始離散點(diǎn)集。首先將幾何域邊界領(lǐng)域（領(lǐng)域半徑δ）的柵格點(diǎn)全部移動(dòng)到邊界上，這樣得到邊界離散點(diǎn)；再刪除幾何域外部的全部柵格點(diǎn)。此時(shí)邊界上的離散點(diǎn)及幾何域內(nèi)的柵格點(diǎn)構(gòu)成了初始離散點(diǎn)集合。

在不更改區(qū)域邊界上節(jié)點(diǎn)位置的前提下，為了提高網(wǎng)格質(zhì)量，應(yīng)對臨近邊界處的內(nèi)部柵格節(jié)點(diǎn)進(jìn)行判斷及優(yōu)化。若幾何域內(nèi)柵格點(diǎn)與邊界的距離大于柵格控制尺寸h，則不需要處理；若柵格點(diǎn)與邊界的距離小于柵格尺寸的一半，則刪除該柵格節(jié)點(diǎn)；若柵格點(diǎn)與邊界的距離d大于柵格尺寸的一半但小于柵格尺寸，則將其向幾何域內(nèi)部移動(dòng)到距離邊界（d＋h）／2處。處理后如圖1所示（圖中h＝0．1，δ＝h／2）。

圖1 簡單幾何域的初始點(diǎn)集

1．2 節(jié)點(diǎn)優(yōu)化

由圖1中可知，幾何域邊界處的網(wǎng)格質(zhì)量差，需要進(jìn)一步的優(yōu)化。利用MATLAB 2013a中的函數(shù)DT＝delaunayTriangulation（p）進(jìn)行Delaunay三角劃分，其中P為n×2階點(diǎn)集坐標(biāo)矩陣，DT為n×3階矩陣，保存構(gòu)成三角單元的三個(gè)頂點(diǎn)的編號(hào)，即P矩陣的行數(shù)。

本文中使用非線性彈性系數(shù)的彈簧比擬法來優(yōu)化網(wǎng)格質(zhì)量。彈簧比擬法的基本思想是將網(wǎng)格節(jié)點(diǎn)與節(jié)點(diǎn)之間的連線看作彈簧，通過彈簧變形產(chǎn)生彈性力移動(dòng)節(jié)點(diǎn)位置，從而實(shí)現(xiàn)對網(wǎng)格質(zhì)量的優(yōu)化。下面說明非線性彈性系數(shù)在節(jié)點(diǎn)位置優(yōu)化中的使用。

由DT可以得到各節(jié)點(diǎn)上的邊，各邊都視為彈簧，故有彈性形變函數(shù)fij（l，l0）＝k（l0－l）（l，l0分別表示實(shí)際長度和所需長度，下標(biāo)ij表示該邊的端點(diǎn)即網(wǎng)格節(jié)點(diǎn)的編號(hào)）來表示各邊受力與邊長的關(guān)系，當(dāng)l＜l0時(shí)，fij＞0，即力的方向沿軸向向外使邊延長；反之，則向內(nèi)使邊縮短。

考慮到網(wǎng)格單元的質(zhì)量，即三角形的三邊長度近似，要求彈性系數(shù)足夠大；再考慮網(wǎng)格劃分的光順，要求當(dāng)實(shí)際長度在指定長度的某一鄰域內(nèi)，即l∈（l0－δ，l0＋δ）時(shí)，彈性系數(shù)應(yīng)足夠小以使節(jié)點(diǎn)便于移動(dòng)，有利于加快節(jié)點(diǎn)位置的優(yōu)化速度。為此引入非線性彈性系數(shù)：

式中，α為加權(quán)系數(shù)；lij和l0分別表示節(jié)點(diǎn)上的任意一邊邊長和所需邊長。利用冪函數(shù)y＝（αx）β的特性，｜x｜＞1／α?xí)r函數(shù)值快速上升；｜x｜＜1／α?xí)r函數(shù)值快速逼近0，尤其是在｜x｜＜1／2α區(qū)域內(nèi)。本文中將曲線斜率變化的膝點(diǎn)定義如下，若取β＝3，要求邊長變化的容許范圍為±20%，則取α≈2．027。如此可以使三角形邊長差之比在區(qū)域（－0．2，0．2）內(nèi)時(shí)剛度系數(shù)較小，使網(wǎng)格尺寸平滑過度；有利于網(wǎng)格頂點(diǎn)達(dá)到受力平衡，減小三角網(wǎng)格各邊之差，在保證網(wǎng)格形狀飽滿、提高均勻度的同時(shí)使網(wǎng)格大小過度順滑。

幾何域邊界頂點(diǎn)不能移動(dòng)，同時(shí)邊界上的其他節(jié)點(diǎn)不能離開邊界，只能沿邊界移動(dòng)。故在幾何域邊界上施加外力，使得邊界結(jié)點(diǎn)不能向外移動(dòng)。

先建立N＊2維矩陣P來存儲(chǔ)有效節(jié)點(diǎn)p的坐標(biāo)。再將各邊的彈力分解為水平分量、垂直分量，故有

Fint表示節(jié)點(diǎn)所受各邊彈力的合力；Fb是人為添加的法向外力，方向向內(nèi)，僅施加在邊界結(jié)點(diǎn)上，以保證邊界結(jié)點(diǎn)只能沿邊界線移動(dòng)。故F（p）的第一、二列分別表示節(jié)點(diǎn)受力的x軸分量和y軸分量。

若各節(jié)點(diǎn)受力平衡，即平衡方程F（p）＝0有解。引入時(shí)間依賴關(guān)系，轉(zhuǎn)變?yōu)槌Ｎ⒎址匠探M［4］，初值為p（0）＝p0，可得

當(dāng)p達(dá)到穩(wěn)態(tài)時(shí)，便滿足方程組F（p）＝0。式（3）可以通過向前Euler方法來近似求解。將時(shí)間離散為tn＝nΔt，可由下式求得近似解pn≈p（tn）

對全部柵格點(diǎn)位置通過迭代計(jì)算優(yōu)化后，再使用Delaunay三角網(wǎng)格劃分法剖分幾何域，所得結(jié)果如圖2所示，邊界處的網(wǎng)格質(zhì)量有明顯提高。

圖2 優(yōu)化后的均勻三角網(wǎng)格剖分

1．3 網(wǎng)格尺寸的控制

單一網(wǎng)格劃分中所需網(wǎng)格尺寸即三角單元的邊長l0為一常量，但很多情況下，幾何域的不同部分對網(wǎng)格大小有不同要求。為此引入單元尺寸函數(shù)h（x，y），給出幾何域中三角單元尺寸的相對系數(shù)。由式（5）可得網(wǎng)格尺寸的變化系數(shù)

式中，h0表示網(wǎng)格最小尺寸，pij為邊lij的中點(diǎn)。f（x，y）可為平面上的點(diǎn)、線或閉合區(qū)域。式（5）表示通過節(jié)點(diǎn)到f（x，y）的距離來控制網(wǎng)格的尺寸。

如圖3所示，左圖中靠近原點(diǎn)處網(wǎng)格尺寸趨近0．025，靠近過點(diǎn)（0．3，1．0）和（1．0，0．5）的直線10x＋14y＋17＝0處的網(wǎng)格尺寸趨近于0．025；右圖中網(wǎng)格尺寸在邊界處趨近0．03。

圖3 不均等網(wǎng)格劃分

2 網(wǎng)格質(zhì)量分析與評估

二維單元的幾何形狀主要是三角形和四邊形，主要質(zhì)量指標(biāo)包括：單元長度、翹曲角、單元邊長比、內(nèi)角大小、扭曲角、雅可比比率（Jacobian ratio）等。三角形單元主要檢查：單元長度、長寬比、扭曲角和內(nèi)角大?。?］。

本文使用常用三角形的形狀因子，即三角單元最大內(nèi)接圓半徑的兩倍與最小外接圓半徑之比［6］

式中，a、b、c表示三角形的邊長；0≤q≤1。若為正三角形，則q＝1；若為退化三角形（即三角形三頂點(diǎn)共線），則q＝0。

通常，若q＞0．5，就認(rèn)為網(wǎng)格質(zhì)量較好，本文中的改進(jìn)方法通常在0．7以上。通常使用Laplacian順滑算法的Delaunay網(wǎng)格劃分的低質(zhì)網(wǎng)格少于4%，而一般的Delaunay限定網(wǎng)格劃分中則高達(dá)10%～20%。本文中所使用的基于非線性彈簧比擬法的網(wǎng)格單元中的高質(zhì)量網(wǎng)格比例更高。如圖4所示，上下兩圖分別顯示圖3中的矩形與圓形的網(wǎng)格單元的質(zhì)量分布。

3 結(jié)束語

圖4 網(wǎng)格質(zhì)量評估

本文方法與一般平面三角網(wǎng)格生成方法相比，其優(yōu)點(diǎn)在于對網(wǎng)格單元各邊引入彈力比擬，以受力平衡來修正網(wǎng)格的質(zhì)量。同時(shí)通過非線性彈性系數(shù)保證網(wǎng)格的光順過度；其缺點(diǎn)在于增加了算法實(shí)現(xiàn)的復(fù)雜性和計(jì)算量，對于大規(guī)模網(wǎng)格十分不利。

［1］ 王成恩．面向科學(xué)計(jì)算的網(wǎng)格劃分與可視化技術(shù)［M］．北京：科學(xué)出版社，2011．

［2］ 古成中，吳新躍．有限元網(wǎng)格劃分及發(fā)展趨勢［J］．計(jì)算機(jī)科學(xué)與探索，2008，2（3）：248－259．

［3］ 江 平，劉民士．射線法判斷點(diǎn)與包含簡單曲線多邊形關(guān)系的完善［J］．測繪科學(xué)，2009，34（5）：220－222．

［4］ Persson P O ，Strang G ．A Simple Mesh Generator in MATLAB［J］．SIAM Review，June 2004，46（2）：329－345．

［5］ 李海峰，吳冀川，劉建波，等．有限元網(wǎng)格剖分與網(wǎng)格質(zhì)量判定指標(biāo)［J］．中國機(jī)械工程，2012，23（3）：368－377．

［6］ 佘紅偉．二維區(qū)域網(wǎng)格剖分算法研究［D］．西安：西北工業(yè)大學(xué)理學(xué)院，2003．