適用于電磁場有限元計算的網格剖分算法

章春鋒 汪 偉 吳天緯 安斯光

(中國計量大學機電工程學院 浙江 杭州 310016)

0 引 言

優異的適應場域邊界幾何形狀以及媒介物理性質變異的能力，使有限元法成為各類電磁場、電磁波工程問題定量分析和優化設計的主導數值計算方法之一[1]。一般而言，有限元分析劃分為建立模型 (前處理)、計算求解、結果處理和評定(后處理)。其中各階段所用的時間占比分別為：40%～60%、5%～10%、30%～50%。因此，如何將實際工程中復雜場域離散為“有限單元”，這一剖分問題在數值分析過程是極其關鍵的一步，直接決定了最終結果的準確性、有效性[2]。

網格生成器諸如大型商業軟件CAE/Maxwell，以映射法為理論基礎[3]。其主要思想是，通過合適的映射函數將待剖分的幾何域映射到參數空間形成規則的參數域；對規則參數域進行網格剖分；將參數域的網格反向映射回幾何域，從而得到幾何域的有限元網格。映射法是非全自動方法，必須通過人工交互方式，將剖分對象先剖分成具有簡單拓撲關系的子域[4]。以映射法為核心算法的這些商業軟件價格昂貴、操作復雜，且在生成大規模三角網格時速度較慢，往往是幾乎不可訪問的復雜代碼，它們通常被用作“黑盒”，很難與其他代碼或者軟件進行交互，因此用戶放棄了控制[5]。理解和使用網格生成算法的能力，好比掌握數據可視化、計算機圖形中的幾何建模的方法一樣，是非常有價值的選擇，值得深入研究[6-7]。

Persson-Strang算法是基于 MATLAB 的有限元網格自動生成算法[8]。與目前常用的Gmsh和Triangle開源網格剖分技術相比，Persson-Strang算法對一般的二維和三維的建模都能進行良好的處理，具有程序簡短明晰、網格質量高、可移植性好的突出優點[9]。但是Persson-Strang算法并不完善，它可能不會終止并返回一個形狀良好的網格，剖分速度相對較慢，且不易與實際應用相結合。目前，網格剖分存在剖分過程耗時長和網格質量不高的問題[10]。

本文主要從四個方面構造新的網格剖分算法：① 通過引入新的平滑函數減少迭代次數；② 利用坐標映射技術，避免由于微小形變而需重新剖分引起的計算資源浪費；③ 加入質量評估來作為終止條件，避免過度迭代；④ 對可能出現的重復節點進行刪除處理。對于有限元應用解決了一些問題：確保邊界節點位于邊界上，提取有限元邊界信息，將三角形節點逆時針編號處理。

1 網格剖分原理

1.1 Persson-Strang算法原理

Persson-Strang算法的簡單性是由于使用有符號距離函數來指定和描述要網格化的幾何區域。距離函數指定從空間中的任何點到域的邊界的最短距離，其中函數的符號在區域外是正的，在內部是負的，在邊界上是零。此定義用于判斷節點是位于幾何區域R2的內部還是位于幾何區域的外部[8]。此外，距離函數d(x)的梯度指向邊界的方向，允許外部的點有效地移回到域。幾何域Ω可以定義為：

Ω={p=(x,y)∈R2|d(p)<0}

(1)

在該算法中，所需的網格大小由網格密度函數h(p)提供。網格密度函數h(p)是由用戶指定的，其大致是區域的邊界之間的距離。在幾何域比較復雜的位置，可以使用自適應方法來實現局部特征變化[11]。對于均勻的網格，h(p)=1。

Persson-Strang算法采用基于力平衡的網格平滑函數來優化節點位置，該力平衡模型是三角形網格和二維桁架結構之間的類比。三角形的邊對應于桿，頂點對應于桁架的節點。平衡問題對應一個非線性方程：

F(p)=0

(2)

找到一個定常解，滿足式(2)。系統使用正向歐拉方法進行近似。在離散時間tk=kΔt時，近似解pk≈p(tk)更新為：

pk+1=pk+ΔtF(pk)

(3)

式中：Δt為方便離散所設置的時間分段值。

(4)

式中：eij=[pi,pj]為可變彈簧常數；f(pi,pj)取決于其當前長度‖eij‖和所需的長度rij。每一條邊的力F取決于它當前的長度‖eij‖和理想長度rij之間的差值。對于均勻網格，rij是常數。但在許多情況下，在不同的區域有不同的尺寸是有利的[12]。輕微的非線性力函數可能會產生更好的網格，但分段線性力給出良好的結果，設置k=1，當‖eij‖=rij時，F=0，這是合理的。建議的邊界處理方法，意味著沒有任何點被迫停留在邊界上，只是被阻止越過邊界。

對于某些幾何形狀(如圓盤、橢球等)，創建帶符號距離函數很容易。對于其他像凸多邊形一樣簡單的幾何域，距離函數很難得到表示。Persson-Strang算法通過引入固定點，可以使用簡化的距離函數。除了比其他網格生成器更短、更簡單之外，Persson-Strang算法還可以生成高質量的網格。

1.2 Persson-Strang可改進點

Persson-Strang算法是一種用于生成非結構化三角形網格的快速而靈活的工具，但是對于有限元應用而言存在一些主要缺點：

(1) 基于Laplacian函數的平滑算法僅描述了排斥力，并未考慮吸引力，造成迭代次數過多，剖分速度緩慢。

(2) 對于有限元方法(FEM)的數值模擬，生成合格的網格既復雜又耗時。在進行實際應用時，即使只有很小的變化，也需要為每組參數重復生成相應的有限元網格。重新網格化過程非常耗時。

(3) 算法執行緩慢和存在不終止的可能性。

(4) 有限元應用方面的不足。不能保證邊界節點位于域的邊界處，固定點和邊界節點在網格剖分排序時可能出現重復。網格剖分核心函數Delaunay剖分生成的網格未逆時針排列編號。

2 改進的網格剖分算法

對于Persson-Strang網格剖分算法的不足之處，本文算法主要給出以下改進。

2.1 平滑過程

為確保有限元分析結果令用戶滿意，生成的有限元網格通常應具備以下條件：所有單元接近理想形狀；變化梯度較大的地方網格密度應相應增大；粗細網格之間過渡均勻[10]。但通常情況下，在有限元網格自動生成器所生成的網格中總存在一些發生畸變的單元，剖分的幾何域越復雜，畸形網格所占比例越大[13]。

本文研究通過平滑函數進行網格優化處理。平滑處理常用的Laplacian函數是應用最早同時也是最成熟的一種優化方法。其核心內容為：保持網格拓撲關系不變，將整個內部節點的位置攝動到由其相鄰節點組成的多邊形的質心處，使每個單元更接近于理想形狀。將這個攝動過程遍歷所有內部節點若干次，可較大幅度地提高網格質量[14]。經過Laplacian光滑處理的網格，不再具有Delaunay三角剖分的基本性質。

此外，化學中的Lennard-Jones函數考慮了吸引與排斥兩種力兩個方面的因素，給出分子間作用能表達式，被稱為蘭納-瓊斯勢。蘭納-瓊斯勢函數是表示分子間相互作用勢能的一種近似模型，如式(5)所示。

f(d)=d-13-d-14

(5)

式中：d是網格元素的邊的長度。

Laplacian平滑函數只描述吸引行為，以幫助將節點分布到整個Ω域，但是由于只具有排斥力，不可避免出現迭代頻繁；蘭納-瓊斯勢函數平滑函數效果受到數值不穩定的影響，有概率導致迭代不終止，陷入死循環[11]。綜合兩個函數的平滑特性，通過擬合方法，從數學上得到新的平滑函數NSF(New Smoothig Function)，NSF同時具有排斥力和吸引力，且在數值上具有穩定性：

f(d)=(1-d2)×ed3-d4

(6)

平滑函數如圖1所示，其中：實點線為NSF平滑算法作用效果；虛線為Laplacian算法作用效果。

圖1 平滑函數

(7)

2.2 自適應微小形變技術

提出一種穩定且簡單的網格變形技術，用于快速準確地劃分FEM網格。首先將邊界網格定義為可伸縮變形的骨架，進而控制網格剖分生成的節點，不包括固定節點。更新所選邊界網格的新坐標，就可以快速映射回細網格中節點的新位置。

在調用Delaunay三角形網格剖分后，進行網格優化之前，通過距離函數nfd(new function distance)定義新的邊界網格框架，通過計算距離函數的數值梯度grad，將點投射到新的邊界。

di(xi,yi)=nfd(xi,yi)

(8)

(9)

式中：deps為預設誤差值，一般為10-3h0。

(10)

(11)

進行坐標映射之后，開始網格平滑優化過程，計算并存儲每個節點的坐標,可通過重心公式判斷是否需要細分，經驗判斷方式為：若重心在三角形內部，則不需要細分，若在外部，則需要細分。該算法既不需要網格重新剖分，也不需要添加新的未知參數。與現有方法相比，計算復雜度及其計算成本較低。該算法使用時，需關注網格質量q值，在形變較大時，仍須重新剖分以確保網格剖分的質量。

2.3 終止條件

該算法存在執行緩慢和不終止的可能性。實驗表明通過附加的控制邏輯，可以使它更加健壯。終止標準基于當前迭代中節點的最大移動小于收斂精度dptol：

(12)

式中：Fr(pi)是每個節點pi的位移量；h0是初始網格尺寸大小。

由于一個高質量的網格有可能在達到該終止條件之前早已形成終止，標準應該包括一個質量評估，以避免迭代頻繁。

在有限元應用中，誤差上限取決于網格中的最小角度。要量化網格質量，通常使用的質量度量是最大內切圓半徑rin與最小外接圓半徑rout之間的關系，即：

(13)

式中：a、b和c是邊長。根據經驗法則，如果所有三角形的qmin>0.5，則網格是優的。在算法中加入該質量評估，以避免過度迭代。

2.4 問題節點處理

指定的固定點與應用概率拒絕方法后幸存下來的點可能重復，這將使得節點編號排序不正確，進而影響有限元求解。

對該問題進行如下處理：移除區域外的點，應用概率拒絕方法，再添加固定節點。在最后網格生成完成時，基于坐標判斷是否重復，一般固定點位于邊界處，在提取邊界節點為有限元做準備后，在有限元解決計算過程中再次檢查是否有重復節點。上述步驟在進行Delaunay三角剖分和平滑優化之前完成。

將三角形的任意兩頂點交換，生成具有負面積的三角形排序，生成的三角形按標準逆時針形式排列；邊界邊在剖分域中只會出現一次，通過判斷邊的出現次數，實現邊界節點篩選。

2.5 算法流程

本文算法流程如下：

步驟1初始化和參數設置：根據初始網格的大小h0，先把能涵蓋欲劃分區域的最大矩形劃分為結構網格。根據距離函數定義，移除邊界外的點。應用概率拒絕方法篩選點，當指定某些點要保留時，將保留的點加入，并刪除重復的點。

步驟2Delaunay檢索：執行Delaunay三角剖分并移除質心在外部的三角形。逆時針重新排列三角形頂點，形成不重復的邊界邊并提取。

步驟3微小形變處理：給出新的距離函數來定義新的邊界框架，計算新距離函數的數值梯度，通過坐標映射技術將點投影到新的邊界上。計算并存儲坐標節點，防止出現節點重復。計算重心，判斷是否需要細分。

步驟4網格平滑：計算和組裝邊緣力并移動節點。

步驟5邊界節點：在邊界上移動網格邊界節點。網格的邊界節點為邊界邊緣的端點。

步驟6終止標準：計算內部點的相對位移，并通過點來確定邊長以計算網格質量q。如果在當前三角剖分中檢測到大位移，則轉到步驟2。如果在當前迭代中檢測到大位移，則轉到步驟3，如果q值均大于0.5，則終止。

步驟7有限元準備：清理并檢查最終網格，輸出包含邊界節點序號、迭代次數等統計信息。

3 驗證分析與應用

3.1 驗證分析

實驗條件：計算機配置為Intel(R)Core(TM)i5-7500 CPU @ 3.40 GHz。

取經典的二維復合幾何網格生成實驗提出的新的平滑函數對于網格剖分迭代次數以及網格生成質量的影響。采用新的平滑函數的網格剖分圖如圖2所示。

圖2 NACA0012剖分

測試Laplacian函數和NSF平滑函數對于Persson-Strang算法的影響，在使用式(3)的前提下，取四組h0進行測試對比，h0為初始定義的網格尺寸大小，結果如表1所示。使用NSF平滑函數優化網格到達平穩狀態所需要的迭代次數是使用Laplacian平滑函數進行網格優化所需要迭代次數的50%。

表1 平滑函數對于迭代次數以及時間影響

通過FEATool Multiphysics實現并集成了統一的MATLAB網格生成框架，基于該工具箱，可以直接比較Persson-Strang算法和本文算法網格生成時間和網格質量(h0=0.02)。

如表2所示，經過四次剖分測試(L1，L2，L3，L4)，新的網格算法相較于原先速度得到提升，其生成的網格質量q值經檢測均在0.7以上。

表2 CPU時間(復合幾何) s

通過二維示例驗證所提出的微小形變處理方法。如圖3所示。假定內部圓形在變形過程中被放大。關于所提出的網格更新方法，首先需要定義包括所有可移動節點的邊界網格。重新定義距離函數，可以根據新的邊界網格(如圖4(a)中黑色虛線所標)和不變的區域坐標來映射精細網格的新位置，得到圖4(b)。

圖3 微小形變示意圖

圖4 確立新的距離函數之后的剖分圖

得到CPU運行時間，以顯示本文方法在計算效率上的優點。此二維精細網格中有399個節點和680個元素。再生網格的CPU時間為7.263 s，而針對相同問題的建議方法的計算時間為6.024 s,縮短了大約17%。兩個網格中每個元素的平均最大與最小邊緣比分別為1.21和1.33。

3.2 應用于電磁場領域

同軸線是由兩根同軸的圓柱導體構成的導行系統，內外導體之間填充空氣或高頻介質的一種寬頻帶微波傳輸線。同軸線主要以TEM模的方式廣泛應用于寬頻帶饋線和元器件的設計中。如圖5所示，使用本文算法剖分一個外導體邊長為2 cm，內導體邊長為1 cm，內外導體電位差為1 V的方同軸線，進行電磁場有限元處理，計算方同軸線間的電位分布，處理步驟如下[15-16]：

圖5 方同軸線圖

(1) 網格剖分。

(2) 讀入網格剖分準備的數據文件。

(3) 計算總體系數矩陣K。

(4) 處理強加邊界條件。

(5) 求解有限元方程。

使用本文算法對該方同軸線進行剖分，在幾何域比較復雜的地方，通過網格密度函數來進行有針對性的細化，網格剖分圖如圖6所示。與Persson-Strang算法的速度和網格質量進行對比，結果如表3所示。

圖6 網格剖分示意圖

表3 不同算法剖分得到的結果

可以觀察到，引入NSF平滑函數后的本文算法，迭代次數是原先Laplacian平滑算法的60.9%左右，運行速度得到顯著提升，得到的網格質量是平穩的，網格質量優，而且相較于Persson-Strang算法，不會出現網格質量顯著下降的情況。

進行電磁場有限元處理后得到的等位線分布圖和電位三維曲面圖如圖7和圖8所示。坐標軸設置單位長度為1 cm，得到以下四點坐標：P1(-0.85,-0.85)，P2(-0.75,-0.75)，P3(-0.65,-0.65)，P4(-0.55,-0.55)。選取四點來驗證電位，如表4所示。

圖7 等位線分布圖

圖8 電位分布三維曲面圖

表4 不同算法對電位求解的影響V

可以看出，本文算法在二維規則形狀的剖分處理上，相較于Persson-Strang算法，速度得到顯著提升，并且維持了優異的網格質量；與Maxwell專業有限元軟件剖分結果相比，本文算法經過有限元后處理返回得到的電磁場參數與商業軟件得到的結果精確度相似，這驗證了新的網格剖分算法網格剖分速度以及質量方面改進的有效性。

4 結 語

本文提出一種新的網格剖分算法，基于Persson-Strang算法網格剖分原理，通過引入新的平滑函數來減少網格迭代優化過程中迭代次數；通過在Delaunay三角形剖分處理后新定義網格邊界框，利用梯度坐標映射技術，避免由于微小形變而須重新進行網格剖分引起的計算機資源浪費；通過設置質量評估來避免網格過度迭代，以此提升網格剖分速度；對網格節點進行編號，并提取邊界信息為有限元處理做準備；最后應用于電磁場領域，驗證了該方法的可行性。