劉 瑤，譚建國

(國防科技大學，湖南 長沙 410073)

0 引 言

一般地，可以將有限元分析[1]劃分為：建立模型(前處理)、計算求解、結果處理和評定(后處理)。其中，各階段所用的時間占比分別為：40%～60%、5%～10%、30%～50%。在數(shù)值分析過程中，網(wǎng)格劃分[2]是極其關鍵的一步，可以說直接決定了最終結果的準確性、有效性。MATLAB[3]具有簡潔的語言、豐富的基礎庫函數(shù)、強大的數(shù)據(jù)處理和可視化能力等突出特點，利用MATLAB進行有限元前處理是值得研究的問題。

DistMesh[4]是基于MATLAB的有限元網(wǎng)格自動生成程序，由加州大學伯克利分校的Per-Olof Persson和麻省理工學院數(shù)學系的Gilbert Strang開發(fā)，相較于MATLAB自帶的PDE工具箱僅能求解比較簡單的二維模型的缺陷，它可以生成三角形或四面體的網(wǎng)格，對一般的二維和三維的建模都能進行良好的處理，具有程序簡短明晰、網(wǎng)格質量高、可移植性好等突出優(yōu)點，在MATLAB有限元前處理方面應用廣泛。

但是由于DistMesh用來篩選和優(yōu)化節(jié)點的距離函數(shù)是解析表達式，對簡單的模型來說可以得到良好的應用，用戶卻無法實現(xiàn)對復雜模型的構造，這限制了DistMesh的工程應用和推廣。目前比較通用的三維建模方式[5]有網(wǎng)格建模、多邊形建模、面片建模、非均勻有理B樣條(NURBS)[6]建模等，而因為良好的數(shù)學特性和強大的建模功能，NURBS能夠比傳統(tǒng)的建模方式更好地控制物體表面的曲線度，從而創(chuàng)建出更逼真、生動的造型，成為CAD中幾何定義的標準技術。通用的三維CAD軟件也都采用NURBS建模方法。

文中基于MATLAB的NURBS工具箱[7]，對DistMesh進行改進，提出一種利用NUBRS曲線和曲面的方向來判斷節(jié)點與曲線、曲面位置關系的方法，以完成對節(jié)點的篩選；再基于力的平衡原理對網(wǎng)格邊的長度進行調整，以完成對節(jié)點位置的迭代優(yōu)化。

1 網(wǎng)格生成器DistMesh

1.1 DistMesh設計原理

DistMesh基本原理是使用距離函數(shù)d(x,y)來定義幾何邊界，通過判斷節(jié)點位置來篩選節(jié)點，再基于力的平衡來完成對網(wǎng)格的迭代優(yōu)化。使用d(x,y)計算每個節(jié)點到邊界的距離，規(guī)定在區(qū)域內部時的距離為負，外部為正，從而將在邊界外的點去除。后續(xù)優(yōu)化是基于三角形網(wǎng)格的邊與桁架結構或者是線性彈簧之間的物理類比[8]。平面中的任何一組點都可以由Delaunay算法[9]進行三角化。在物理模型中，三角形網(wǎng)格的邊(頂點之間的連接)對應于桁條，頂點對應于桁架的節(jié)點。

每一條邊的力F取決于它的當前長度l及其理想長度l0的差值：

(1)

網(wǎng)格的疏密分布與幾何形狀的復雜程度有關[10]。在該算法中，所需的網(wǎng)格大小由網(wǎng)格密度函數(shù)h(x,y)提供。h(x,y)由用戶指定，可以使用自適應方法[11]來實現(xiàn)局部特征變化，其大致是區(qū)域的邊界之間的距離。自適應求解器可估計每個三角形中的誤差，從而給出約束條件來精細化網(wǎng)格以獲得更好的分布。

1.2 二維網(wǎng)格生成器算法描述

function[p,t]=distmesh2d(fd，fh，h0，bbox，pfix，varargin)

程序的輸入?yún)?shù)如下：

(1)幾何形狀由距離函數(shù)fd給出。d=fd(p)返回值為每個節(jié)點位置P到其最近邊界的帶符號距離。

(2)所需的邊緣長度h(x,y)由函數(shù)fh給出，為所有輸入點返回h。

(3)h0是網(wǎng)格的初始大小。對于均勻網(wǎng)格h(x,y)=C。

(4)邊界(最大外圍矩形范圍)。

(5)固定節(jié)點(必須出現(xiàn)在網(wǎng)格中的點)。

輸出有：

(1)節(jié)點位置P：該N×2數(shù)組包含N個節(jié)點的x,y坐標。

(2)三角指標t：指定三角形的節(jié)點編號。

基本步驟為：

(1)在涵蓋待劃分幾何形狀的邊界框中劃分初始分布，初始分布可由用戶任意規(guī)定。

(2)利用距離函數(shù)除去不在邊界內的所有節(jié)點。

(3)進入主循環(huán)，節(jié)點位置不斷迭代優(yōu)化：

①去除質心在邊界外的三角網(wǎng)格；

②計算理想的三角形邊長，基于與實際邊長的差值調整邊長、移動節(jié)點；

③將移動到邊界外的點拉回到邊界上，避免節(jié)點流失。

(4)終止循環(huán)，終止條件由當前迭代中的節(jié)點最大移動距離決定。

1.3 存在的問題

距離函數(shù)是解析表達的，對于復雜圖形很難構造，這大大限制了DistMesh的適用性。概括來說，算法中需要改進的地方有：步驟2，利用距離函數(shù)去除邊界外的節(jié)點；進入循環(huán)后的第1步，用Delaunay方法劃分為網(wǎng)格后，去除質心在邊界外的三角網(wǎng)格；進入循環(huán)后的第3步，將移出的節(jié)點移回邊界。

2 基于NURBS的網(wǎng)格生成方法

2.1 改進方法

(2)

圖1 點的位置判斷

同理，對于三維實體，如圖2所示的球面，其法線方向是u方向和v方向的向量積。曲面法線朝外，從曲面外一點pout到其最近點pntm引向量，由于最近點pntm處法線和所引向量的夾角大于90°，二者數(shù)量積必為負；反之，對于曲面內的點數(shù)量積應為正。同理，當存在內表面時，若內表面法線方向朝內，則內表面內的點(即邊界外的點)向最近點所引向量與最近點處法線的數(shù)量積為負，邊界內的點為正。故可以定義外表面法線方向朝外，內表面法線方向朝內。

圖2 通過點與曲面位置關系判斷其內外

2.2 算法描述

二維平面指定調用函數(shù)：function[p,t]=nrbdistmesh2d(crvout,crvin,fh,h0,bbox,pfix)，函數(shù)的輸入?yún)?shù)由fd變成了crvout和crvin；用戶只需要使用NUBRS工具箱進行內外環(huán)曲線的建模。

三維實體指定調用函數(shù)：function[p,t]=nrbdistmesh3d(srfout,srfin,fh,h,box,fix),同理，用戶只需要使用NUBRS工具箱進行內外曲面的建模。

基本步驟如下：

(1)先對涵蓋用戶待劃分區(qū)域的外圍矩形(長方體)創(chuàng)建初始分布：二維平面為三角形網(wǎng)格，三維實體為四面體網(wǎng)格。

(2)判斷外環(huán)曲線(曲面)的方向：基于上述原理，先在邊界外確定一點，通過該點在邊界外這一既定結論反推曲線(曲面)的方向。

對于三維實體，從外圍長方體表面上找一點pout，這一點必在邊界外；找到外表面上最近點pntmout，求出pntmout所在處法線方向wout；再從pout向pntmout引向量，計算這個向量與法線方向的數(shù)量積zout；若zout小于0，說明外表面法線方向朝外，反之，改變外表面方向。

(3)判斷用戶是否輸入內環(huán)(內表面)：如果存在，應當先定義其方向；否則直接執(zhí)行第5步。

(4)定義內環(huán)曲線(內表面)的方向：內環(huán)曲線為順時針方向，內表面法向向內。

對于平面，與外環(huán)原理相同，將外圍矩形上的點替換為剛剛得到的外環(huán)上的A點，此點必在內環(huán)外；當這個點指向最近點的向量與最近點處導數(shù)的向量積z分量為正，說明內環(huán)是順時針；否則，變換方向。對于三維實體，也從外表面上找點，符合要求的內表面，得到的數(shù)量積應該是正的，否則將內表面變換方向。

(5)除去邊界外的點。與原方法計算距離的正負不同，這里仍然利用向量積(數(shù)量積)。

對于平面，因為滿足條件的點都在曲線左側或在曲線上，故而向量積的Z分量小于0的都是在邊界外的，需要去除；對于三維實體，因為外表面法線朝外，內表面法線朝里，滿足條件的點都在曲面法線方向的相反側，則從節(jié)點向邊界上最近點引向量，與最近點處的法線方向是相同的(小于90°)，故而數(shù)量積小于0的都是需要去除的。

(6)進入循環(huán)，進行節(jié)點的迭代優(yōu)化。

①利用網(wǎng)格質心與內外邊界的位置關系，去除質心在邊界外的網(wǎng)格；

②計算理想的網(wǎng)格邊長，基于與實際邊長的差值調整邊長、移動節(jié)點；

③將移動到邊界外的點拉回到邊界上，把邊界外點相對應的邊界上最近點賦給邊界外點即可。

(7)終止迭代循環(huán)，與原方法相同。

3 驗證與分析

3.1 網(wǎng)格質量評價指標

對于三角形網(wǎng)格，網(wǎng)格質量[12]常根據(jù)三角形最大內切圓與其最小外接圓的半徑比率來判斷：

(3)

其中，a,b,c是三角形的邊長。對于正三角形來說，q=1；對于退化三角形，q=0。如果所有的三角形網(wǎng)格都有q>0.5，可認為是高質量的。

3.2 二維平面改進結果

(1)單位圓和帶孔正五邊形的優(yōu)化過程如圖3和圖4所示，從左至右分別是迭代1次、10次、100次的效果。從圖中可以看出，邊界處的三角形網(wǎng)格變化尤為明顯，邊界也趨于清晰理想化，最終整體網(wǎng)格質量很高，q都在0.99以上。

圖3 單位圓網(wǎng)格優(yōu)化過程

圖4 帶孔正五邊形網(wǎng)格優(yōu)化過程

3.3 三維實體改進結果

圖5為球體的四面體網(wǎng)格優(yōu)化過程的剖面圖，由形狀千奇百怪最后趨于均勻的正四面體。

圖5 球體四面體網(wǎng)格優(yōu)化過程

3.4 結果分析

以上所有二維的例子q的平均值超過0.98，顯著好于拉普拉斯平滑的Delaunay細化網(wǎng)格[14]。和原來的DistMesh進行對比，如表1所示，其網(wǎng)格質量仍然保持在高水平。

表1 改進算法與DistMesh網(wǎng)格質量對比

4 結束語

文中介紹了基于MATLAB的網(wǎng)格自動生成器DistMesh的主要原理，由于其核心的距離函數(shù)是解析表達的，對于復雜圖形很難構造，故基于MATLAB中的NUBRS工具箱，將NURBS建模嵌入到算法中，解決了DistMesh不能實現(xiàn)復雜模型構造的缺陷。在二維平面的改進中，利用點到邊界曲線上最近點處的向量，與最近點處切線方向向量的向量積，來定義曲線方向、判斷點與曲線的位置關系，以篩選節(jié)點；在三維實體的改進中，利用點到邊界曲面上最近點處的向量，與最近點處法線方向向量的數(shù)量積，來定義曲面的法線方向、判斷點與曲線的位置關系，以篩選節(jié)點；基于力的平衡原理的迭代循環(huán)，通過對比理想的與實際的網(wǎng)格長度來不斷調整網(wǎng)格邊長，完成對網(wǎng)格的優(yōu)化。驗證結果表明，改進算法在保證網(wǎng)格生成質量的前提下，提高了DistMesh的適用性。