數(shù)值流形法局部網(wǎng)格加密算法

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

(海南大學(xué)土木建筑工程學(xué)院， ?？?570228)

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

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

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

數(shù)值流形法使用局部網(wǎng)格加密以提高精度的方法穩(wěn)定高效，近年來研究局部網(wǎng)格加密方法不是特別多。現(xiàn)通過對避免出現(xiàn)懸掛節(jié)點的加密方案進行研究并利用C++面向?qū)ο缶幊碳夹g(shù)開發(fā)程序，完善了數(shù)值流形法的前處理方案，提高數(shù)值流形法的計算精度。

1 數(shù)值流形覆蓋系統(tǒng)

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

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

圖1 數(shù)值流形法的覆蓋系統(tǒng)Fig.1 Cover system of numerical manifold method

圖2 數(shù)學(xué)覆蓋Fig.2 Mathematical cover

圖3 數(shù)學(xué)覆蓋生成物理覆蓋Fig.3 Generation of physical cover from mathematical cover

圖4 數(shù)學(xué)覆蓋中的流形單元Fig.4 Manifold elements in a typical mathematical cover

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

圖6 一個流形單元關(guān)聯(lián)的物理覆蓋Fig.6 Physical Covers associated with a manifold element

以上三個概念之間的相互關(guān)系對程序設(shè)計以及進行局部加密算法非常重要。本文相關(guān)聯(lián)的程序開發(fā)以Visual Studio為平臺，使用C++面向?qū)ο蟪绦蛟O(shè)計，數(shù)學(xué)覆蓋，物理覆蓋以及流形元分別使用MathCover、 PhysicalCover和ManifoldElement三個類進行定義。程序設(shè)計中使用了C++STL中Map和List等數(shù)據(jù)結(jié)構(gòu)，采用指針變量保存數(shù)學(xué)覆蓋、物理覆蓋以及流形單元。

2 局部加密方案

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

2.1 一到四級局部網(wǎng)格加密

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

圖7 數(shù)學(xué)網(wǎng)格局部加密示意圖Fig.7 Schematic diagram of local refinement of mathematical grid

2.2 局部加密算法

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

圖8 局部網(wǎng)格加密流程圖Fig.8 Local grid refinement flow chart

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

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

(3)根據(jù)程序預(yù)設(shè)的加密半徑，若數(shù)學(xué)網(wǎng)格中心點至裂紋尖端的距離在加密半徑內(nèi)，則標(biāo)記該數(shù)學(xué)網(wǎng)格為需被加密的網(wǎng)格，若數(shù)學(xué)網(wǎng)格中心點至裂紋尖端距離超出加密半徑之外，該數(shù)學(xué)網(wǎng)格無需被加密。

(4)根據(jù)加密等級，需要逐一對被加密的數(shù)學(xué)網(wǎng)格添加加密點及加密線，對需加密區(qū)域進行加密。

(5)為防后續(xù)生成流形單元中含面積過小單元，調(diào)整距離節(jié)理線過近的數(shù)學(xué)網(wǎng)格點到其在節(jié)理線的垂足點上。

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

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

2.3 裂紋尖端應(yīng)力強度因子

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

與路徑無關(guān)的J積分[16],可以表示為

(1)

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

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

(2)

(3)

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

3 算例驗證

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

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

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

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

該問題的應(yīng)力強度因子SIF解析解[20]為

(4)

(5)

劃分網(wǎng)格生成的流形單元如圖10所示，此時網(wǎng)格尚未加密，用已經(jīng)編制完成的二維數(shù)值流形計算程序進行計算分析，具體計算方法可參考計算矩形板在y方向的位移uy云圖并由Tecplot軟件繪制[1-2]，如圖11所示。在裂紋尖端分別進行一至四級局部加密，裂紋尖端加密區(qū)域如圖12所示，紅色虛線框內(nèi)為加密區(qū)域；局部加密后求解得到的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 加密裂紋尖端加密區(qū)域Fig.12 Crack tip refinement area of refinement

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

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

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

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

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

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

該問題的應(yīng)力強度因子SIF解析解為

(6)

(7)

利用已經(jīng)設(shè)置好的計算程序求解該問題，劃分網(wǎng)格生成的流形單元如圖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

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

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

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

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

4 結(jié)論

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

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

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

(3)數(shù)學(xué)網(wǎng)格加密方案分為一至四級，等級越高網(wǎng)格加密越精細(xì)。符合加密越精細(xì)，計算結(jié)果越精確的規(guī)律，說明本文算法有效。