非結構網格自適應細化的實現與應用?

李于鋒

（中國工程物理研究院 計算機應用研究所，四川 綿陽 621900）

0 引言

在工程數值模擬領域，針對復雜結構的有限元計算應用十分廣泛。在機械零件制造中如何高效得到結構關鍵點的應力，以及在流場中如何快速準確地捕捉激波等問題，若單純依靠不斷全局細化的網格增加計算量，是不經濟的。網格自適應技術兼顧了計算精確性和計算效率，能夠使用較少的網格計算代價獲得較為準確的計算結果［1］。

自適應網格技術在國外發展十分迅速，在成熟的商業軟件中都有功能體現，如LS－Dyna的重劃分網格的自適應實現等；在美國Sandia實驗室開發的多物理耦合計算框架SIERRA［2］中，也實現了H自適應的策略。在國內，北京應用物理與計算數學研究所開發的JASMIN［3］結構化網格自適應框架已經在多個應用程序中得到實際應用。中科院科學與工程計算國家重點實驗室開發的PHG［4］能夠對四面體網格進行自適應細化。然而涵蓋二維和三維非結構化網格自適應的軟件模塊還相對缺乏，本文采取基于拼片修復的誤差估計方法，嘗試在自主研發的有限元計算框架中實現二維三角形網格和三維四面體網格的自適應細化，在多層細化中考慮網格質量，并用算例驗證該方法的有效性。

1 后驗誤差估計

后驗誤差估計是通過對有限元計算結果進行再處理以提高計算精度的后驗方法，這種過程稱為修復。修復解更接近精確解，以此作為標準衡量有限元解的誤差的方法稱為后驗誤差估計。誤差定義為精確解和近似解的差，對于位移u，定義為：

其中：L為線性微分算子，p為已知函數。

定義能量范數為：

其中：D為彈性矩陣；S為應變關于位移的算子。由算子S定義的應變和應力分別為：

將式（1）代入式（4），聯立式（5）、式（6），能量范數也可以寫為：

誤差能量范數中的真實應力是未知的，采用超收斂的拼片修復（Superconvergent Patch Recovery）［5，6］方法獲得較準確的應力或應變，來對有限元近似解進行衡量。

2 三角形網格細化策略及實現

為了保持網格質量，二維三角形網格的自適應細化采取互連三個邊中點將單元一分為四的正則細化方案。在細化單元和非細化單元之間會產生過渡單元，這些單元的一條邊或者兩條邊被細化。一個單元在細化過程中可能被標識細化邊的情形如圖1所示。

圖1 三角形可能的細化情況

圖1（a）為正則細化情形，將三條邊的中點互連即可將原單元一分為四；圖1（c）的單元只有一條邊被細化，直接將細化邊中點與相對的頂點連接，可將單元一分為二；圖1（b）的單元有兩條邊被細化，此類單元細化要考慮質量問題，細化方式有如圖2所示的兩種選擇。經過簡單計算，就可以選擇較好的細分方式。

圖2 過渡單元的質量控制

三角形網格經過多次細化后，一個單元可能會產生多次非正則細化情形，從而導致較差的網格質量，為了保證多層細化之后，網格質量不至于太差，需要對細化過程進行一定的約束控制。為此設定一個規則：若將被細化單元本身不是正則細化而得到的，則返回到該單元的父單元，進行正則細化后，再細化該單元（可能是正則或非正則細化）。該準則控制網格質量在實際中取得較好效果，如圖3所示。

圖3 三角形網格細化的質量控制

3 四面體網格細化策略及實現

三維四面體網格由于對復雜幾何模型更容易逼近，在工程實際的建模中應用廣泛。然而四面體網格的細化由于維度的增加，比三角形的細化算法更加復雜。首先，若選擇每條邊細化的完全細化方式，可以將單元一分為八（見圖4），四個頂點處形成4個四面體外，內部還可剖分形成4個四面體，然而內部的小四面體并不能保證質量，由于網格協調性的約束，導致過渡網格的剖分情形復雜而且質量更差。為此，本文采取二分法作為四面體的細分方案。二分法實現簡單，一次細化只需添加一個新節點，也不至于產生過多的單元，見圖5。

圖4 四面體的八分細化 圖5 四面體的二分細化

為了在程序中實現二分細化，根據網格幾何協調性總結出了三種細化方法，分別是邊細化法、逐單元細化法和細化邊列表法。邊細化法是被標細化單元全部細化最長邊，若某單元有多個邊被細化的，則按照邊的長度進行順序二分。該方法的缺點是一次細化將導致某單元多次細化，網格質量很差。逐單元細化法是被標識細化單元即時細化最長邊，但若由于鄰居細化已經細化了，則此次就不再繼續細化它的子樹，該方法導致短邊細化的情形較多且質量不好。細化邊列表法是將所有被標識待細化單元的最長邊存入一優先隊列，按照長度遞減排序，遍歷該細化邊列表，對每個細化邊，即時細化該邊所在的單元，為了不至于產生過多的細化，將引起該邊被細化的單元記為該細化邊的主單元，若遍歷細化邊列表時，該邊的主單元（可能有多個）都已經被細化，則跳過該細化邊。經改進后的細化邊列表法能在實際算例中表現出網格質量的極大改善。

4 自適應網格實例分析

4.1 L形結構板的應力分析

經典的L形結構的物理模型如圖6（a）所示，頂部和右側為滑移邊界，左側受均布拉力，材料為彈性。該結構在L的拐角處容易產生應力集中，自適應細化的網格很好地展示了該問題，用較少的網格量達到了較為精確的解，相對能量范數誤差從初始網格的24%下降到最終自適應網格的4%，初始網格和自適應網格分別見圖6（b）和圖6（c）。

4.2 帶孔板的拉伸分析

無限大的平面板受單向拉伸的物理模型如圖7（a）所示，利用對稱性和圣維南原理，模擬計算其四分之一的原始網格見圖7（b），經過兩次自適應細化之后的網格見圖7（c），相對能量范數誤差從原始的10%降到4%。

4.3 三維懸臂短梁的應力分析

懸臂短梁受頂端均布壓力的物理模型如圖8（a）所示，原始四面體網格見圖8（b），采用質量控制后的自適應網格見圖8（c），相對能量范數誤差從45%下降 到7%。

圖6 L形結構的應力分析

圖7 平板拉伸分析

圖8 三維懸臂短梁應力分析

5 結論

在大尺度復雜結構的有限元模擬中，非結構網格自適應是一種關鍵技術，它能夠在保證精度的前提下，大大縮減計算代價。本文算例中的自適應網格細化能僅在增加局部細化工作量的情況下，獲得較好的模擬效果，表明了非結構網格自適應算法的正確性和有效性。

［1］Zienkiewicz O C，Taylor R L，Zhu J Z.The finite element method：Its basis and fundamentals ［M］.6th ed.England：Elsevier Ltd，2005：456－524.

［2］Stewart J R，Edwards H C.A framework approach for developing parallel adaptive multiphysics applications［J］.Finite Elements in Analysis and Design，2004，40：1599－1617.

［3］Mo Zeyao，Zhang Aiqing，Cao Xiaolin，et al.JASMIN：aparallel software infrastructure for scientific computing［J］.Front Comput Sci China，2010，4（4）：480－488.

［4］Zhang Linbo.PHG：A toolbox for developing parallel adaptive finite element programs ［J］.Bulletin of the Chinese Academy of Sciences，2011，25（4）：298－300.

［5］Zienkiewicz O C，Zhu J Z.Superconvergent patch recovery and a posteriori error estimation in the finite element method，Part 1：A general superconvergent recovery technique［J］.Internat J Num Meth Eng，1992，33：1331－1364.

［6］Zienkiewicz O C，Zhu J Z.The superconvergent patch recovery（SPR）and a posteriori error estimates，Part 2：Error estimates and adaptivity［J］.Internat J Num Meth Eng，1992，33：1365－1382.