李于鋒

（中國工程物理研究院 計(jì)算機(jī)應(yīng)用研究所，四川 綿陽 621900）

0 引言

在工程數(shù)值模擬領(lǐng)域，針對復(fù)雜結(jié)構(gòu)的有限元計(jì)算應(yīng)用十分廣泛。在機(jī)械零件制造中如何高效得到結(jié)構(gòu)關(guān)鍵點(diǎn)的應(yīng)力，以及在流場中如何快速準(zhǔn)確地捕捉激波等問題，若單純依靠不斷全局細(xì)化的網(wǎng)格增加計(jì)算量，是不經(jīng)濟(jì)的。網(wǎng)格自適應(yīng)技術(shù)兼顧了計(jì)算精確性和計(jì)算效率，能夠使用較少的網(wǎng)格計(jì)算代價(jià)獲得較為準(zhǔn)確的計(jì)算結(jié)果［1］。

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

1 后驗(yàn)誤差估計(jì)

后驗(yàn)誤差估計(jì)是通過對有限元計(jì)算結(jié)果進(jìn)行再處理以提高計(jì)算精度的后驗(yàn)方法，這種過程稱為修復(fù)。修復(fù)解更接近精確解，以此作為標(biāo)準(zhǔn)衡量有限元解的誤差的方法稱為后驗(yàn)誤差估計(jì)。誤差定義為精確解和近似解的差，對于位移u，定義為：

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

定義能量范數(shù)為：

其中：D為彈性矩陣；S為應(yīng)變關(guān)于位移的算子。由算子S定義的應(yīng)變和應(yīng)力分別為：

將式（1）代入式（4），聯(lián)立式（5）、式（6），能量范數(shù)也可以寫為：

誤差能量范數(shù)中的真實(shí)應(yīng)力是未知的，采用超收斂的拼片修復(fù)（Superconvergent Patch Recovery）［5，6］方法獲得較準(zhǔn)確的應(yīng)力或應(yīng)變，來對有限元近似解進(jìn)行衡量。

2 三角形網(wǎng)格細(xì)化策略及實(shí)現(xiàn)

為了保持網(wǎng)格質(zhì)量，二維三角形網(wǎng)格的自適應(yīng)細(xì)化采取互連三個(gè)邊中點(diǎn)將單元一分為四的正則細(xì)化方案。在細(xì)化單元和非細(xì)化單元之間會產(chǎn)生過渡單元，這些單元的一條邊或者兩條邊被細(xì)化。一個(gè)單元在細(xì)化過程中可能被標(biāo)識細(xì)化邊的情形如圖1所示。

圖1 三角形可能的細(xì)化情況

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

圖2 過渡單元的質(zhì)量控制

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

圖3 三角形網(wǎng)格細(xì)化的質(zhì)量控制

3 四面體網(wǎng)格細(xì)化策略及實(shí)現(xiàn)

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

圖4 四面體的八分細(xì)化 圖5 四面體的二分細(xì)化

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

4 自適應(yīng)網(wǎng)格實(shí)例分析

4.1 L形結(jié)構(gòu)板的應(yīng)力分析

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

4.2 帶孔板的拉伸分析

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

4.3 三維懸臂短梁的應(yīng)力分析

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

圖6 L形結(jié)構(gòu)的應(yīng)力分析

圖7 平板拉伸分析

圖8 三維懸臂短梁應(yīng)力分析

5 結(jié)論

在大尺度復(fù)雜結(jié)構(gòu)的有限元模擬中，非結(jié)構(gòu)網(wǎng)格自適應(yīng)是一種關(guān)鍵技術(shù)，它能夠在保證精度的前提下，大大縮減計(jì)算代價(jià)。本文算例中的自適應(yīng)網(wǎng)格細(xì)化能僅在增加局部細(xì)化工作量的情況下，獲得較好的模擬效果，表明了非結(jié)構(gòu)網(wǎng)格自適應(yīng)算法的正確性和有效性。

［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.