基于非線性彈簧近似的平面三角網格優化方法

2014-01-01 02:54:42張思全

通信電源技術 2014年5期

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

（上海海事大學物流工程學院，上海201306）

0 引言

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

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

由于Delaunay三角劃分算法可以使全部三角形的最小內角為最大，且保證三角單元不會重疊（即所有邊最多只被兩個三角形共有），同時三角單元的三邊之上不會有其他節點，故三角剖分以Delaunay法的應用最為廣泛［2］。

1 算法思想與實現步驟

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

對于不夠飽滿的單元形狀通過彈簧比擬對節點進行受力分析的迭代計算來使節點達到受力平衡，從而提高單元網格的質量。將網格節點視為受力點Pi（i＝1，2，3，．．．，n），每個Pi都有對應的力－位移函數。受力不平衡（連接到節點的邊的長度不等）使節點產生相應的位移，同時利用MATLAB中的Delaunay三角形函數調整拓撲結構。最終邊長應接近所要求的相對長度（相對長度為1時，則各邊長度幾乎相等）。

1．1 初始離散點集生成

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

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

圖1 簡單幾何域的初始點集

1．2 節點優化

由圖1中可知，幾何域邊界處的網格質量差，需要進一步的優化。利用MATLAB 2013a中的函數DT＝delaunayTriangulation（p）進行Delaunay三角劃分，其中P為n×2階點集坐標矩陣，DT為n×3階矩陣，保存構成三角單元的三個頂點的編號，即P矩陣的行數。

本文中使用非線性彈性系數的彈簧比擬法來優化網格質量。彈簧比擬法的基本思想是將網格節點與節點之間的連線看作彈簧，通過彈簧變形產生彈性力移動節點位置，從而實現對網格質量的優化。下面說明非線性彈性系數在節點位置優化中的使用。

由DT可以得到各節點上的邊，各邊都視為彈簧，故有彈性形變函數fij（l，l0）＝k（l0－l）（l，l0分別表示實際長度和所需長度，下標ij表示該邊的端點即網格節點的編號）來表示各邊受力與邊長的關系，當l＜l0時，fij＞0，即力的方向沿軸向向外使邊延長；反之，則向內使邊縮短。

考慮到網格單元的質量，即三角形的三邊長度近似，要求彈性系數足夠大；再考慮網格劃分的光順，要求當實際長度在指定長度的某一鄰域內，即l∈（l0－δ，l0＋δ）時，彈性系數應足夠小以使節點便于移動，有利于加快節點位置的優化速度。為此引入非線性彈性系數：

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

幾何域邊界頂點不能移動，同時邊界上的其他節點不能離開邊界，只能沿邊界移動。故在幾何域邊界上施加外力，使得邊界結點不能向外移動。

先建立N＊2維矩陣P來存儲有效節點p的坐標。再將各邊的彈力分解為水平分量、垂直分量，故有

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

若各節點受力平衡，即平衡方程F（p）＝0有解。引入時間依賴關系，轉變為常微分方程組［4］，初值為p（0）＝p0，可得

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

對全部柵格點位置通過迭代計算優化后，再使用Delaunay三角網格劃分法剖分幾何域，所得結果如圖2所示，邊界處的網格質量有明顯提高。

圖2 優化后的均勻三角網格剖分

1．3 網格尺寸的控制

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

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

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

圖3 不均等網格劃分

2 網格質量分析與評估

二維單元的幾何形狀主要是三角形和四邊形，主要質量指標包括：單元長度、翹曲角、單元邊長比、內角大小、扭曲角、雅可比比率（Jacobian ratio）等。三角形單元主要檢查：單元長度、長寬比、扭曲角和內角大小［5］。

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

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

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