章春鋒 汪 偉 吳天緯 安斯光

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

0 引 言

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

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

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

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

1 網(wǎng)格剖分原理

1.1 Persson-Strang算法原理

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

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

(1)

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

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

F(p)=0

(2)

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

pk+1=pk+ΔtF(pk)

(3)

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

(4)

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

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

1.2 Persson-Strang可改進(jìn)點

Persson-Strang算法是一種用于生成非結(jié)構(gòu)化三角形網(wǎng)格的快速而靈活的工具，但是對于有限元應(yīng)用而言存在一些主要缺點：

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

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

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

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

2 改進(jìn)的網(wǎng)格剖分算法

對于Persson-Strang網(wǎng)格剖分算法的不足之處，本文算法主要給出以下改進(jìn)。

2.1 平滑過程

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

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

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

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

(5)

式中：d是網(wǎng)格元素的邊的長度。

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

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

(6)

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

圖1 平滑函數(shù)

(7)

2.2 自適應(yīng)微小形變技術(shù)

提出一種穩(wěn)定且簡單的網(wǎng)格變形技術(shù)，用于快速準(zhǔn)確地劃分FEM網(wǎng)格。首先將邊界網(wǎng)格定義為可伸縮變形的骨架，進(jìn)而控制網(wǎng)格剖分生成的節(jié)點，不包括固定節(jié)點。更新所選邊界網(wǎng)格的新坐標(biāo)，就可以快速映射回細(xì)網(wǎng)格中節(jié)點的新位置。

在調(diào)用Delaunay三角形網(wǎng)格剖分后，進(jìn)行網(wǎng)格優(yōu)化之前，通過距離函數(shù)nfd(new function distance)定義新的邊界網(wǎng)格框架，通過計算距離函數(shù)的數(shù)值梯度grad，將點投射到新的邊界。

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

(8)

(9)

式中：deps為預(yù)設(shè)誤差值，一般為10-3h0。

(10)

(11)

進(jìn)行坐標(biāo)映射之后，開始網(wǎng)格平滑優(yōu)化過程，計算并存儲每個節(jié)點的坐標(biāo),可通過重心公式判斷是否需要細(xì)分，經(jīng)驗判斷方式為：若重心在三角形內(nèi)部，則不需要細(xì)分，若在外部，則需要細(xì)分。該算法既不需要網(wǎng)格重新剖分，也不需要添加新的未知參數(shù)。與現(xiàn)有方法相比，計算復(fù)雜度及其計算成本較低。該算法使用時，需關(guān)注網(wǎng)格質(zhì)量q值，在形變較大時，仍須重新剖分以確保網(wǎng)格剖分的質(zhì)量。

2.3 終止條件

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

(12)

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

由于一個高質(zhì)量的網(wǎng)格有可能在達(dá)到該終止條件之前早已形成終止，標(biāo)準(zhǔn)應(yīng)該包括一個質(zhì)量評估，以避免迭代頻繁。

在有限元應(yīng)用中，誤差上限取決于網(wǎng)格中的最小角度。要量化網(wǎng)格質(zhì)量，通常使用的質(zhì)量度量是最大內(nèi)切圓半徑rin與最小外接圓半徑rout之間的關(guān)系，即：

(13)

式中：a、b和c是邊長。根據(jù)經(jīng)驗法則，如果所有三角形的qmin>0.5，則網(wǎng)格是優(yōu)的。在算法中加入該質(zhì)量評估，以避免過度迭代。

2.4 問題節(jié)點處理

指定的固定點與應(yīng)用概率拒絕方法后幸存下來的點可能重復(fù)，這將使得節(jié)點編號排序不正確，進(jìn)而影響有限元求解。

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

將三角形的任意兩頂點交換，生成具有負(fù)面積的三角形排序，生成的三角形按標(biāo)準(zhǔn)逆時針形式排列；邊界邊在剖分域中只會出現(xiàn)一次，通過判斷邊的出現(xiàn)次數(shù)，實現(xiàn)邊界節(jié)點篩選。

2.5 算法流程

本文算法流程如下：

步驟1初始化和參數(shù)設(shè)置：根據(jù)初始網(wǎng)格的大小h0，先把能涵蓋欲劃分區(qū)域的最大矩形劃分為結(jié)構(gòu)網(wǎng)格。根據(jù)距離函數(shù)定義，移除邊界外的點。應(yīng)用概率拒絕方法篩選點，當(dāng)指定某些點要保留時，將保留的點加入，并刪除重復(fù)的點。

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

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

步驟4網(wǎng)格平滑：計算和組裝邊緣力并移動節(jié)點。

步驟5邊界節(jié)點：在邊界上移動網(wǎng)格邊界節(jié)點。網(wǎng)格的邊界節(jié)點為邊界邊緣的端點。

步驟6終止標(biāo)準(zhǔn)：計算內(nèi)部點的相對位移，并通過點來確定邊長以計算網(wǎng)格質(zhì)量q。如果在當(dāng)前三角剖分中檢測到大位移，則轉(zhuǎn)到步驟2。如果在當(dāng)前迭代中檢測到大位移，則轉(zhuǎn)到步驟3，如果q值均大于0.5，則終止。

步驟7有限元準(zhǔn)備：清理并檢查最終網(wǎng)格，輸出包含邊界節(jié)點序號、迭代次數(shù)等統(tǒng)計信息。

3 驗證分析與應(yīng)用

3.1 驗證分析

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

取經(jīng)典的二維復(fù)合幾何網(wǎng)格生成實驗提出的新的平滑函數(shù)對于網(wǎng)格剖分迭代次數(shù)以及網(wǎng)格生成質(zhì)量的影響。采用新的平滑函數(shù)的網(wǎng)格剖分圖如圖2所示。

圖2 NACA0012剖分

測試Laplacian函數(shù)和NSF平滑函數(shù)對于Persson-Strang算法的影響，在使用式(3)的前提下，取四組h0進(jìn)行測試對比，h0為初始定義的網(wǎng)格尺寸大小，結(jié)果如表1所示。使用NSF平滑函數(shù)優(yōu)化網(wǎng)格到達(dá)平穩(wěn)狀態(tài)所需要的迭代次數(shù)是使用Laplacian平滑函數(shù)進(jìn)行網(wǎng)格優(yōu)化所需要迭代次數(shù)的50%。

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

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

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

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

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

圖3 微小形變示意圖

圖4 確立新的距離函數(shù)之后的剖分圖

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

3.2 應(yīng)用于電磁場領(lǐng)域

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

圖5 方同軸線圖

(1) 網(wǎng)格剖分。

(2) 讀入網(wǎng)格剖分準(zhǔn)備的數(shù)據(jù)文件。

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

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

(5) 求解有限元方程。

使用本文算法對該方同軸線進(jìn)行剖分，在幾何域比較復(fù)雜的地方，通過網(wǎng)格密度函數(shù)來進(jìn)行有針對性的細(xì)化，網(wǎng)格剖分圖如圖6所示。與Persson-Strang算法的速度和網(wǎng)格質(zhì)量進(jìn)行對比，結(jié)果如表3所示。

圖6 網(wǎng)格剖分示意圖

表3 不同算法剖分得到的結(jié)果

可以觀察到，引入NSF平滑函數(shù)后的本文算法，迭代次數(shù)是原先Laplacian平滑算法的60.9%左右，運行速度得到顯著提升，得到的網(wǎng)格質(zhì)量是平穩(wěn)的，網(wǎng)格質(zhì)量優(yōu)，而且相較于Persson-Strang算法，不會出現(xiàn)網(wǎng)格質(zhì)量顯著下降的情況。

進(jìn)行電磁場有限元處理后得到的等位線分布圖和電位三維曲面圖如圖7和圖8所示。坐標(biāo)軸設(shè)置單位長度為1 cm，得到以下四點坐標(biāo)：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

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

4 結(jié) 語

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