數值流形法局部網格加密算法

伍書重， 莫思陽， 張友良， 張亞軍

(海南大學土木建筑工程學院， 海口 570228)

數值流形法[1-2](numerical manifold method，NMM)通過有限覆蓋技術，使用數值流形覆蓋系統劃分網格，相比有限元方法可以更加精確地計算連續及非連續物體，主要在網格劃分、覆蓋形式、近似函數方面有巨大優勢。

使用數值流形法進行具體模型求解時，一般有兩種方法可以提高計算精度。一種方法是通過使用高階的覆蓋函數來達到提高計算精度的目的。田榮等[3]開發一、二階流形方法程序說明高階覆蓋函數提高計算精度的有效性；蔡永昌等[4-5]使用四節點四邊形的數值流形覆蓋系統并使用高階覆蓋位移函數提高精度；鄧安福等[6]混合使用高、低階覆蓋函數，使得求解的精度和效率均得到提高；郭朝旭等[7]根據近似函數的剛度矩陣提出了LDLT算法[L指下三角單位矩陣(lower triangular identity matrix)，D指對角矩陣(diagonal matrix)，LT指L的轉置矩陣(transpose matrix)]，穩定快速求得特解。高階近似函數帶來的高精度算法會出現線性相關的問題[8]。

另一種提高計算精度的方法就是網格加密。很多學者采用整體網格加密方法[9-11]，然而整體加密方法雖然提高了精度，但是卻嚴重降低了計算效率，增加了計算時間。具體問題分析時，只有部分區域需要進行網格加密，張友良等[12]提出基于等幾何分析的數值流形法，基于T樣條局部加密數學覆蓋網格，既提高了計算精度，又考慮了計算效率。還有一些學者借鑒有限元成熟的前處理方法進行局部網格加密[13-15]。

數值流形法使用局部網格加密以提高精度的方法穩定高效，近年來研究局部網格加密方法不是特別多。現通過對避免出現懸掛節點的加密方案進行研究并利用C++面向對象編程技術開發程序，完善了數值流形法的前處理方案，提高數值流形法的計算精度。

1 數值流形覆蓋系統

數值流形方法相對于有限元方法具有的優勢之一是數值流形覆蓋系統，數值流形覆蓋系統由三部分組成：數學覆蓋(mathematical cover，MC)、物理覆蓋(physical cover，PC)和流形單元(manifold elements，ME)。數學覆蓋記為MCi(i=1，2，…，nMC)。物理覆蓋由數學覆蓋和物理網格組成，物理覆蓋是數學覆蓋的再剖分，由分析模型的邊界線、裂縫、材料線等組成，物理覆蓋記為Pij。流形元是數學覆蓋和物理覆蓋的交集，流形元記為eij。

數學覆蓋由一個個數學網格構成，可由用戶選擇數學覆蓋的形狀以及疏密程度，但是一般根據問題需要選擇規則的網格。數學網格(如圖1中的3×6四邊形網格)記為mi(i=1,2,…,nm)，數學網格中數學網格線的交點就是數學覆蓋中心點(①～⑨)，如圖2中圍繞數學覆蓋中心點的4個矩形構成一個數學覆蓋MC。

圖1 數值流形法的覆蓋系統Fig.1 Cover system of numerical manifold method

圖2 數學覆蓋Fig.2 Mathematical cover

圖3 數學覆蓋生成物理覆蓋Fig.3 Generation of physical cover from mathematical cover

圖4 數學覆蓋中的流形單元Fig.4 Manifold elements in a typical mathematical cover

圖5 物理覆蓋中的流形單元Fig.5 Manifold element in physical cover

圖6 一個流形單元關聯的物理覆蓋Fig.6 Physical Covers associated with a manifold element

以上三個概念之間的相互關系對程序設計以及進行局部加密算法非常重要。本文相關聯的程序開發以Visual Studio為平臺，使用C++面向對象程序設計，數學覆蓋，物理覆蓋以及流形元分別使用MathCover、 PhysicalCover和ManifoldElement三個類進行定義。程序設計中使用了C++STL中Map和List等數據結構，采用指針變量保存數學覆蓋、物理覆蓋以及流形單元。

2 局部加密方案

局部加密是計算精度與計算效率的協同。網格加密屬于數值流形方法前處理，在例如裂紋尖端這樣需要計算精度高處加密局部數學網格，在原有的數學網格布置完成以后，生成數值流形法覆蓋系統并判斷需要加密的數學網格區域，加密以后，加密區域原來生成的數學覆蓋、物理覆蓋以及流形單元往往會變化，需要更新以保證數值流形覆蓋系統正確。局部更新數學覆蓋、物理覆蓋和流形元，形成新的計算單元并進行計算分析。

2.1 一到四級局部網格加密

面對不同的計算實例，不同的計算模型有不同的精度需求，本研究提出數學網格一級到四級的局部加密算法來滿足不同的精度需求。該算法將數學網格加密等級分為一至四級，數學網格最終一級到四級的加密形態如圖7所示。首先在裂紋尖端周圍生成一個合適大小的四邊形，再將四邊形的4個頂點與對應數學網格的頂點連接形成一級局部加密網格如圖7(a)所示；分別在頂點之間的四條連接線二等分并且將對應的二等分點連接形成二級局部加密網格如圖7(b)所示；同理，分別在頂點之間的四條連接線三等分和四等分并且將對應點連接即可形成三級局部加密網格和四級加密網格，如圖7(c)、圖7(d)所示。

圖7 數學網格局部加密示意圖Fig.7 Schematic diagram of local refinement of mathematical grid

2.2 局部加密算法

對一個生成完成的數值流形覆蓋系統的數學網格加密后，所產生新的網格已經不是適合分析的數值流形覆蓋系統，要保證數值流形法的覆蓋系統適用，必須在原有數學覆蓋系統的基礎上進行數學覆蓋、物理覆蓋以及流形元的更新，這也是本算法的難點所在。局部加密算法流程圖如圖8所示，該算法主要步驟如下：

圖8 局部網格加密流程圖Fig.8 Local grid refinement flow chart

(1)讀取模型數據文件，根據計算精度需要確定加密等級。

(2)識別需要加密的區域，如裂紋尖端、應力集中區等應力梯度較大區域。

(3)根據程序預設的加密半徑，若數學網格中心點至裂紋尖端的距離在加密半徑內，則標記該數學網格為需被加密的網格，若數學網格中心點至裂紋尖端距離超出加密半徑之外，該數學網格無需被加密。

(4)根據加密等級，需要逐一對被加密的數學網格添加加密點及加密線，對需加密區域進行加密。

(5)為防后續生成流形單元中含面積過小單元，調整距離節理線過近的數學網格點到其在節理線的垂足點上。

(6)根據加密后數學網格，依次生成加密后的數學覆蓋MC、流形單元ME和物理覆蓋PC。

(7)由于采用四邊形進行網格劃分，檢查流形單元ME所對應的數學覆蓋MC和物理覆蓋PC數量是否均為4個，若未輸出錯誤信息則加密完成。

2.3 裂紋尖端應力強度因子

利用已經設計好的數值流形方法程序計算裂紋尖端應力強度因子(stress intensity factor，SIF)與理論數值解進行比較是判斷計算精度的有效方法。在斷裂力學理論中，應力強度因子K是預測斷裂發生和裂紋在帶裂紋構件中擴展速率的主要標準[16]。計算應力強度因子的方法有相互作用積分法[17]、J積分法、位移外推法等。本研究采用J積分法計算應力強度因子。

與路徑無關的J積分[16],可以表示為

(1)

式(1)中：Γ為按逆時針方向環繞裂紋尖端的回路曲線，曲線起點為裂紋下表面，終點為裂紋上表面；W為曲線Γ任一點處應變能密度；tx和ty為曲線Γ任一點處的應力分量；ux和uy分別為曲線Γ任一點處的位移分量；ds為沿回路曲線的弧單元。

線彈性條件下，對Ⅰ型純張拉裂紋，J積分與KI之間存在一定的比例關系，平面應力、應變情況可以表示為

(2)

(3)

數值流形方法模擬不連續問題時出現的裂紋尖端奇異性問題，參照劉登學等[18]提出的通過將附加函數加入到覆蓋函數中解決，使用這種方法，裂紋尖端可落在流形單元的任意位置。即使在相對粗糙的數學網格中，也可以精確地得到相應的應力強度因子。

3 算例驗證

3.1 含有中心邊裂紋的有限矩形板拉伸

如圖9所示有限板[19]，側面伸進一條長為a的邊裂紋，有限板單向受拉力σ，矩形板長為2h，寬度為b，其他力學參數如表1所示。

表1 含中心邊開裂紋矩形板的物理力學參數Table 1 Physical and mechanical parameters of rectangular plate with single edge crack

圖9 受單向拉伸的帶單邊裂紋有限板Fig.9 Finite plate with single edge crack in tension

該問題的應力強度因子SIF解析解[20]為

(4)

(5)

劃分網格生成的流形單元如圖10所示，此時網格尚未加密，用已經編制完成的二維數值流形計算程序進行計算分析，具體計算方法可參考計算矩形板在y方向的位移uy云圖并由Tecplot軟件繪制[1-2]，如圖11所示。在裂紋尖端分別進行一至四級局部加密，裂紋尖端加密區域如圖12所示，紅色虛線框內為加密區域；局部加密后求解得到的uy圖如圖13所示。

圖10 含中心國邊裂紋矩形板生成的流形單元Fig.10 Generated manifold elements in rectangular plate with sigle edge crack

圖11 含中心邊裂紋矩形板y方向位移云圖Fig.11 uy contour plots of rectangular plate with sigle edge crack

圖12 加密裂紋尖端加密區域Fig.12 Crack tip refinement area of refinement

代入具體數據后得到SIF解析解為2.877。裂紋尖端局部區域分級加密后，求解得到相應的SIF值，各級加密求得的SIF值及其與解析解的誤差如表2所示。從表2中可以看出，不加密和局部加密后得到的SIF與解析解誤差均不大，說明加密算法求解是有效可行的；同時，從計算結果可以看出局部加密等級越高，誤差越小，得到的SIF值越精確，這也證明了此一到四級加密算法的可行性。

表2 帶單邊裂紋矩形板不同加密條件下的SIFTable 2 SIF of rectangular plate with single edge crack under different refinement conditions

3.2 有限板雙邊裂紋受單向拉伸

圖14所示為一受單向拉伸的帶雙邊裂紋有限板[19]。該矩形板底部邊緣y方向位移被約束，且左下角節點位移被x、y雙向約束。該算例滿足平面應力條件。矩形板的各物理力學參數如表3所示。

圖14 受單向拉伸的帶雙邊裂紋有限板Fig.14 Finite plate with double edge cracks in tension

表3 帶雙邊裂紋矩形板的物理力學參數Table 3 Physical and mechanical parameters of rectangular plate with double edge cracks

該問題的應力強度因子SIF解析解為

(6)

(7)

利用已經設置好的計算程序求解該問題，劃分網格生成的流形單元如圖15所示，Tecplot軟件繪制出矩形板在y方向的位移如圖16所示。

圖15 帶雙邊裂紋矩形板生成的流形單元Fig.15 Generated manifold elements in rectangular plate with double edge cracks

圖16 帶雙邊裂紋矩形板y方向位移云圖Fig.16 uy contour plots of rectangular plate with double edge cracks

在裂紋尖端分別進行一至四級局部加密，各級加密區域及加密網格與3.1算例類似。圖17所示為各級局部加密后求解得到的uy圖。

圖17 加密帶雙邊裂紋矩形板y方向位移云圖Fig.17 uy contour plots of refinement rectangular plate with double edge cracks

代入具體數據后得到SIF解析解為1.199。裂紋尖端局部區域分級加密后，利用已設計好的數值流形法計算程序解得到相應的SIF值，各級加密求得的SIF值及其與解析解的誤差在表3、表4中列出。從表4中可以看出，不加密和局部加密后得到的SIF與解析解誤差均不大，說明利用已開發的計算程序求解是有效可行的；同時，局部加密等級越高，得到的SIF值越精確，這也證明了此加密算法的可行性。

表4 帶雙邊裂紋矩形板不同加密條件下的SIFTable 4 SIF of rectangular plate with double edge cracks under different refinement conditions

4 結論

在非線性數值流形的基礎上，考慮了模擬帶裂紋結構時在裂紋尖端位移場和應力場不精確的問題，引入附加插值函數對奇異物理覆蓋位移函數進行完善，并采用與路徑無關的J積分方法計算應力強度因子SIF，提出了針對加密區域數學網格的加密方案予以解決。得出結論如下：

(1)對兩個經典的帶裂紋平板算例進行了計算分析，按照不加密、一級加密、二級加密、三級加密、四級加密的順序對裂紋尖端需加密區域進行網格加密，程序計算得到對應的應力強度因子SIF，并列出不同加密等級下求得的SIF與解析解之間的誤差，兩算例存在的誤差均較小，各級加密方案所得SIF誤差均在5%以內。

(2)隨著加密等級的提高，SIF與解析解間存在誤差趨于降低，在相同半徑的J積分圓軌跡條件下，第四級加密誤差最小。

(3)數學網格加密方案分為一至四級，等級越高網格加密越精細。符合加密越精細，計算結果越精確的規律，說明本文算法有效。