基于Delaunary三角剖分的WOA優化算法研究

周 燁，吳 芳，溫 瑋，李沛宗，楊日杰

(海軍航空大學，山東 煙臺 264100)

實現水下目標聲散射精確建模是檢驗反潛聲吶檢測算法的重要環節。針對水下目標進行有限元數值仿真可以精確模擬水下目標聲散射現象，從而有效仿真水下聲場。在有限元計算中，大約60%的工作在于前期處理。其中，有限元網絡劃分是非常關鍵的，網絡剖分的質量將直接影響最終的計算結果準確性。目前針對二維平面剖分，Delaunary三角剖分方法由于具有較好的理論基礎和數學屬性，已經成為有限元分析軟件前處理的廣泛選擇。Delaunary剖分算法的主要優勢體現在空圓特性以及最大化最小角理論，其避免了狹長三角形單元的產生[1]。可以極大地避免低質量單元格的產生。但是在實際應用中仍有很多問題，比如在復雜狹長區域的劃分不均，域外三角單元問題等。近年來，針對Delaunary剖分，有很多學者提出了優化方案，比如L.A.Piegl提出的邊界外法向法，凸分解法等，在解決域外三角單元方面都有各自的優缺點[2]。針對這種情況，文獻[1]提出了一種算法優化方案，較好地解決了域外三角單元和劃分不均的問題。但是在處理實際問題中，特別是艦艇水聲問題下的有限元仿真建模，我們往往不要求求解域中的單元均勻劃分，在流固耦合作用區域，我們希望加密網絡得到較為精確的解，在水體等大范圍環境參數變化不大的區域，我們希望減少單元格來提高運算速度。為了解決這個問題，本文中利用權重賦值建立艦艇艙室結構不同重點區域，在流固耦合界面、固體力學傳遞界面和聲波傳遞界面建立不同的加密區域，利用控制單元節點插入次序在不同區域建立均等三角單元網絡，提高網絡整體質量[3]。

1 基本原理

1.1 經典Delaunary三角算法

先對Delaunary算法進行簡單介紹。針對其方法的實現途徑有很多，其中B-W算法因操作簡單、易于實現等優點被普遍應用[4]。其核心思想是循環利用插入點與已知點性質來生成細化三角形[5]。其具體操作步驟如下：

步驟1根據輸入模型進行識別，生成初始離散點；

步驟2在邊界點的基礎上構造包圍劃分區域的初始三角單元；

步驟3根據B-W算法依次構建已知離散點與初始三角單元相互作用的小三角單元；

步驟4刪除包含有初始三角單元的三角形，得到已知點的三角剖分；

步驟5根據優化準則進行局部優化；

步驟6循環執行步驟3。其中B-W算法步驟如圖1所示。

圖1 B-W剖分步驟示意圖

1.2 艦艇艙室結構的權重比例因子

在水聲仿真主動探測領域，艦艇復雜的艙室結構特征會給回波信號帶來不可忽略的影響，相互鏈接的肋狀支撐結構從一定程度上使艦艇艙室的本征振動頻率發生改變，從而出現意料之外的共振峰[6]，精確模擬艙室結構特征對于聲吶回波信號仿真及目標識別具有重要意義。本文以潛艇結構為例，其結構模型也可以對艦船模型構建產生良好的借鑒。某型雙殼潛艇構建內部艙室結構模型如圖2所示，結構細化如圖3所示。

圖2 某型潛艇艙室結構模型示意圖

圖3 結構細化示意圖

由于部分區域涉及流固耦合以及聲固相互作用，需要進行網格加密來提高計算精度，[7]部分區域因為單元數值變化梯度低，可以采用稀疏網格單元增加計算速度，因此需要對不同網絡作用區域進行權重調整。首先，以重要程度劃分，對水體和大型艙室空氣作用權重賦初始值[8]，其余結構依次是肋板接觸區域、單一流固接觸界面、聲固耦合作用再輻射區域。通過全加密剖分下的有限元應力分析，可以大致指定權重比例。采用模型整體最下網絡剖分下，潛艇艙室結構橫切面模型對1 000 Hz平面波水平激勵下的結構應力分析如圖4所示。依據此可以大體指定艙室內部的權重值配比。

針對應力作用區域，劃分權重值，如表1所示。

圖4 應力分析示意圖

表1 權重劃分

1.3 區域優化算法

由于Delaunary剖分算法在局部區域進行細化插值時，插值點選取不同而導致低質量網格出現，如圖5所示。

圖5 隨機點選取示意圖

當△ABC和△ABE為待細分三角單元時，傳統算法可能會在△ABC中插入點，這樣會造成圖6中的情況，產生了△ADC與△BDE兩種單元，這是有限元分析中不想見到的。

圖6 傳統算法生成網絡示意圖

在相同類型區域計算時，希望盡可能用一組等邊三角形去擬合網絡。[9]因此，在相同權重區域進行剖分時，可以在待細分單元中選擇可以生成最佳三角單元的區域插入新節點，并實時檢驗有效性。其算法實施思路總結如下：

步驟1識別模型邊界，并進行均勻離散，生成初始的Delaunary三角單元集合p1。

步驟2判斷域外三角形。計算所有邊界三角形的重心D，設其三角形的邊界邊分別為ci和cj，按照右手定則判斷域外三角形，若方向為正向，則刪除該單元。

步驟3針對權重比例，設置不同的網格加密區域并設置相關期望網格尺寸{αi}，其中i表示權重等級，αi表示不同等級下的期望尺寸。在權重區域判斷上，可以重點區域的邊界點所構造的初始三角形為區域細化目標，針對不同區域設置不同集合，下面的細化操作均在不同區域的集合內進行操作。

步驟4首先針對最低等級進行全模型網格細化。通過for循環，對所有初始網格進行尺寸對比。設目標三角形邊長最大值為M，對滿足M>α1的所有三角單元進行分層。將(α1,M)區間劃分為小型區間，其過程可以表示為：

(1)

式(1)中：Z表示動態變量；I是隨著每次三角形質量檢查時的Z的遞變量；L表示對單元質量進行檢查的次數。

步驟5將滿足條件M>Z=(Z-I·i)的三角形單元確定為有限插入單元，并添加至統一的集合p2中。

步驟6當一次循環中有多個三角形滿足步驟4中條件時，需要確定插入節點的順序。這里利用邊長排序進行區分。將這些三角形按最長邊降序排列，如果出現最長邊相等，就以邊長之和大的靠前排序；若三邊之和也相等，就比較其最長邊與最短邊的比值，較大的靠前。

步驟7將新節點插入排序后的三角形單元的形心處，構造新的三角單元并刪除原有單元。之后利用Delaunary原理進行優化并更新集合p1和p2。

步驟8當滿足條件的單元節點插入結束后，更新i的值，進行下一階段的for循環，在下個階段的循環中，Z的值也隨之變換。當i>L時，結束這一區域的加密。跳轉至下一區域重復步驟4之后的計算處理。

還是以圖4中三角形為例，在細分算法優化后，點會先在△ABE中產生，從而得到較高質量單元，相對分布較為均勻，如圖7所示。

圖7 均分后網絡

2 算例演示

利用Matlab對Delaunary剖分進行算法優化。先利用Comsol對潛艇初始橫截面模型進行構建，保存成matlab可以直接讀取的dxf文件。進行初始網格劃分，如圖8所示。

圖8 Matlab初始網格

隨后根據設置的閾值進行細化插值，其改善前后艙室細化結構如圖9、圖10所示。

圖9 傳統局部細化

圖10 均分局部細化

由于模擬實際尺度計算，僅給出局部來清晰展示剖分效果。由圖9、圖10可以發現，優化算法對于由Delaunary剖分帶來的網格不均勻現象起到了明顯的改善效果。

3 算法檢驗

3.1 網格質量

針對三角形網絡質量，仿照文獻[1]，引入網絡平均質量系數qm進行判定，即：

(2)

(3)

式(3)中：ri為三角形內切圓半徑；Ri為三角形外接圓半徑，當其趨于1時，為最優情況。對比結果如表2所示。

表2 檢驗對比

結果表明，WOA優化算法可在整體上提高了Delaunary剖分的網格質量[10]。繪制網絡質量直方圖，如圖11所示。

圖11 傳統算法(上)與優化算法(下)網絡質量統計直方圖

由圖11可知，WOA優化算法相對傳統算法無論是總體質量還是質量分布上均有改善，有效滿足了有限元仿真計算的需求。同時發現，優化算法較傳統算法在網格數量上有所增加，原因可能是：為了保證計算精度，在部分區域設置的較低期望值使得優化算法在部分區域產生了比傳統算法更多的細化網絡所致。下面討論兩者在計算精度以及計算時間上的差別。

3.2 計算精度及效率

為了檢驗優化算法剖分后的計算精度，以權重最高的區域期望值為最終期望值，進行全模型極細化剖分計算，以追求最大限度的準確性，并與利用優化算法后針對同一模型的相同位置處的散射聲壓級的仿真計算結果作圖，如圖12所示。

圖12 仿真計算結果

由圖12可以發現，傳統有限元算法在常規密度下的計算結果與極細化的剖分結果在低頻上相似度較高，但是在高頻處有所偏離。原因是沒有考慮不同結構帶來的計算影響，部分重點計算區域沒能細化處理，使得結果發生偏離。而改進后的算法在計算結果上與極細化剖分后的的算法在全頻段保持一致，其網格數以及計算時間如表3所示。

表3 計算效率統計

由表1可以發現，在達到極細化計算精度的情況下，通過建立相應的權重等級優化可以極大減少構建網格數量，降低仿真計算時間。另外，優化算法雖然在網格數量上較傳統算法有一定的劣勢，但是在計算時間上相差不大，綜合計算精度考慮，優化算法在計算復雜結構體方面是可取的。

4 結論

本文提出了針對艦艇復雜艙室結構水聲有限元仿真的WOA算法優化，也可以拓展至其他領域的復雜結構計算。針對實際模型仿真涉及流固耦合以及聲固相互作用，需要進行網格加密提高計算精度的需求設計權重配比，再對相同權重區域進行有效優化剖分，WOA算法可以在保證計算精度的情況下降低畸形網絡生成概率，提高網絡質量和計算速度。本文中的優化算法僅僅涉及二維平面計算，如何將其拓展到三維區域進行有效構建是下一步工作的重點。