等幾何分析與近場動力學耦合的殼結構裂紋擴展仿真算法

夏陽， 王洪帥， 鄭國君， 申國哲

(大連理工大學 工業裝備結構分析國家重點實驗室 汽車工程學院，遼寧 大連 116024)

0 引 言

斷裂仿真對殼結構的安全評估具有重要意義。近年來，斷裂力學的數值仿真取得快速發展，為業界提供多種可行的解決方案，如內聚力單元法、擴展有限元方法、擴展等幾何方法、無網格法、相場方法、分子動力學方法和近場動力學方法等。目前，斷裂仿真存在的問題是如何對結構中連續場和非連續場進行統一描述，并且在計算中保證計算效率。為解決這一問題，斷裂力學非連續模型與連續模型相結合是一個值得研究的方向。

等幾何分析(isogeometric analysis, IGA)直接將精確幾何用于數值分析，連通幾何設計與數值分析之間的鴻溝，適用于殼結構分析的等幾何列式不斷被提出，包括無轉動的Kirchhoff-Love(K-L)殼模型、帶轉動自由度的Reissner-Mindlin(R-M)殼模型、實體殼等。基于裂紋擴展問題，將擴展有限元中用到的單位分解法引入IGA中，得到基于K-L理論和R-M理論的擴展IGA(extended IGA, XIGA)方法。XIGA能夠模擬裂紋附近的應力場，但是需要額外的準則處理裂紋成核和擴展，在處理裂紋分叉和合并方面存在一定困難。另外一種處理斷裂問題的耦合方法是IGA與相場方法耦合。相場方法是建立在Griffith能量斷裂準則之上的，該準則是定量描述斷裂擴展的經典準則，其核心是假設裂紋擴展過程中的總能量(表面能和體積能)保持最小。AMBATI等采用等幾何實體殼單元與相場方法耦合的方式處理殼結構的脆性和韌性斷裂。KIKIS等將考慮剪切變形的等幾何R-M殼與相場方法相結合，避免體積離散化。為能夠得到合適的結果，相場方法必須在裂紋擴展處進行網格細化以滿足收斂的要求，導致其需要消耗大量的計算時間。

近場動力學(peridynamics, PD)理論由SILLING提出，用于處理非連續場問題。學者們提出很多用于板殼結構問題的PD模型。TAYLOR等通過將3D模型降為2D，提出PD薄板模型，該模型能夠捕獲面外彎曲變形，但是只能處理具有特定泊松比的材料。DIYAROGLU等提出考慮剪切變形的薄/厚梁和板的PD列式，其運動方程可以簡化為經典Timoshenko梁和Mindlin板方程。SHEN等基于微梁鍵提出6自由度Euler梁和K-L殼模型，通過插值方法獲得鍵的變形、建立鍵的微勢能，從而得到鍵的微模量，該模型不限制泊松比。以此為基礎，考慮剪切變形，進一步發展出PD R-M殼模型。

由于PD是一種非局部理論，每一個點與其近場域內的一群點相互作用，計算效率比經典的有限元或IGA更低。采用非局部理論和經典連續介質力學理論耦合的方式，其計算效率有望得到提高。將IGA與PD耦合是進行裂紋擴展仿真的一種新的重要方法。MADENCI等提出IGA與PD耦合處理平面問題的方法，將PD點加入到等幾何模型中，PD點的自由度與等幾何節點相關聯，但是PD點不出現在最終的方程中。XIA等提出用于平面彈性問題的IGA-PD耦合模型，使用等幾何控制點作為PD點，根據力平衡方法將IGA方程與PD方程相疊加得到待求解線性方程，是一種簡單且高效的求解方法。此外，鄭國君等使用神經網絡預測裂紋擴展路徑，快速且準確。

殼的光滑曲面特征使其適于采用IGA方法進行模擬。為對殼結構進行裂紋擴展分析，本文提出一種用于板殼結構斷裂分析的耦合模型IGA-PD。采用IGA模型保證幾何精確性，將PD模型用于裂紋擴展區域。為直觀顯示裂紋擴展過程，采用非連續伽遼金(discontinuous Galerkin, DG)方程，提出DG有限元PD板殼單元列式。與粒子點法PD列式相比，采用DG有限元PD板殼單元可以使用高斯點積分，提高計算精度。在耦合過程中，采用力/力矩平衡原理得到耦合模型。將IGA作為邊界可以消除PD的表面效應問題，提高PD模型性能。

1 R-M殼的IGA和微梁鍵模型

1.1 R-M殼IGA模型

采用由三維彈性模型退化得到的R-M殼模型，其幾何模型示意見圖1。殼變形前采用直法線假設，其物理場采用非均勻有理B樣條(non-uniform rational B-spline, NURBS)插值，表達式為

圖 1 R-M殼模型示意

(,,)=(,)+·(,)=

(1)

式中：和為節點參數,可分別為由節點向量=[…++1]和=[…++1]定義；∈[-1,1]為殼的厚度；=[]，即殼上控制點的位置矢量；(,)為與控制點有關的二維張量積NURBS基函數；為控制點處殼體中面的法向量，=[，1，2,3]，采用Greville橫標點計算。

基于等參概念，采用與幾何場相同的插值函數，殼的位移表示為

(,,)=

(2)

對于階基函數，完全積分需要+1個高斯點。位移場與轉動場不一致會導致剪切閉鎖，因此使用個高斯點的減縮積分方案避免這一問題。

通過最小位能原理，最終得到等幾何R-M殼單元的剛度矩陣

=?d=

(3)

在圖1中給定一點,其坐標為(,,)，根據式(2)，點的速度可以通過點計算得到，即

=+

(4)

其中：

(5)

點處體積微元的動能為

(6)

根據直法線假設，點相對于點的位置表示為=[0 0]。將式(4)代入式(6)，得到采用中面表示的動能方程，即

(7)

式中：為材料密度；

(8)

(9)

采用等幾何形函數進行離散，對式(7)在單元域上進行積分，得到一致質量矩陣，即

(10)

1.2 適用于R-M殼的微梁鍵

近場動力學是一種使用積分方程求解連續和非連續場的理論。該理論假設粒子與周圍一定距離范圍內的其他粒子存在相互作用，損傷發生在粒子水平上。參考構型中2個粒子的相對位置為

=′-

(11)

相對位移為

=(′,)-(,)

(12)

這2個相互作用的粒子之間形成PD鍵。粒子作用在粒子上的力可由微勢能求導得到，即

(13)

微勢能是單根鍵上的能量，其單位是單位體積的平方，所以給定點的單位體積的能量即局部應變能密度

(14)

系數12是指每根鍵的每個端點上只有鍵上12的能量。基于近場動力學理論的微梁鍵，將PD鍵視為類似于經典Timoshenko梁模型，承受軸向變形、扭轉變形、彎曲變形和橫向剪切變形，稱之為微梁鍵，見圖2。

圖 2 PD微梁鍵示意

每個微梁鍵的端點有3個平動自由度(,,)和3個轉動自由度(,,)，通過插值方式可獲得鍵的變形。鍵在局部坐標系下的位移向量為

(15)

根據文獻[14-15]，為描述殼結構的變形，微梁鍵上一點的變形可以分成面內變形和橫向變形兩部分，即

,Shell=,PlaneShell+,TranShell

(16)

面內變形微勢能由端點平動產生的軸向變形能和端點轉動產生的橫向變形能組成，則

(17)

式中：a,為鍵在拉伸和壓縮作用下的微模量；b,為與面內彎曲相關的微模量；為鍵的長度。

由殼橫向變形相關的微勢能可得

(18)

式中：為與轉動相關的微模量；b,和s,分別為由于鍵繞軸轉動產生的與彎曲和橫向剪切相關的微模量。

根據應變能密度守恒原理，微梁鍵基殼模型的應變能密度必須與彈性殼理論的應變能密度相等。通過插值離散微梁鍵模型，得到微模量的表達式為

(19)

最終，局部坐標系下微梁鍵的力密度與位移的關系表示為

=

(20)

(21)

(22)

的具體細節可參閱文獻[14-15]。

2 非連續伽遼金PD殼

在經典PD模型中，物理域被離散成粒子，稱為無網格PD離散方法；近幾年發展的另一種方法是將模型域離散成單元，再應用DG方程。例如，REN等使用DG模型離散裂紋存在的非局部模型，使用有限元法離散連續模型域。與粒子法相比，DG方法可以在計算域內使用更高精度的積分，可以更方便地計算應力和應變場。采用微梁鍵基PD殼模型計算的系統總能量為

(23)

(24)

式中：為全局坐標系下點和的位移；為全局坐標系到鍵局部坐標系下的旋轉矩陣。

使用DG方法時，需要先將模型離散成單元，本文采用由控制網形成的4節點單元，見圖3。為所有單元的共享節點插入新的節點，解除單元之間的節點關聯。這種單元的拆分過程可保證各個單元的獨立性，使裂紋擴展分析可以順利進行。

圖 3 DG離散中的單元拆分示意

在PD殼的DG法中，單元是計算模型的基本組成部分。在列式中，單元的高斯點被視為PD節點，見圖4。

圖 4 DG單元的高斯點之間形成PD鍵

單元的高斯點被視為附著在單元上的PD節點。兩單元的任意2個高斯點之間形成鍵，如圖中鍵連接點和點，點的近場域由單元決定。單元中的4個PD點都采用單元形心形成的近場域，單元和單元之間的關聯通過鍵的相互作用建立，即

(25)

令和分別表示單元坐標系與全局坐標系之間位移向量的轉換矩陣，則

(26)

將式(24)～(26)代入式(23)，得

(27)

采用DG法的單元與單元之間的單剛矩陣可寫為

=

(28)

對所有單元的剛度矩陣進行累加，即得到結構的整體剛度矩陣。

3 基于力/力矩平衡原理的耦合模型

DG列式方法增加單元節點的數量，因此降低計算效率。為提高裂紋擴展仿真的計算效率，根據前文獲得的IGA殼模型和DG殼的PD模型，采用非局部模型和經典連續介質力學模型耦合的方法，得到殼結構的動力學平衡方程，即

(29)

方程的每一行表示一個自由度所處的一種力力矩平衡狀態。處于一種邊界條件下的結構，其所承受的等效載荷向量相等，即滿足=。當采用IGA模型或PD模型求解時，2種模型得到的位移解應該相同，即滿足=。因此，方程又可表示為

(30)

根據式(30)，可提出基于力力矩平衡的耦合求解方法。將計算域按等幾何控制點形成的網格劃分為IGA域與PD域，見圖5。

圖 5 耦合模型示意

由此，耦合剛度矩陣可表示為

(-())+()

(31)

式中：為對角矩陣，

(32)

且滿足

(33)

耦合后的動力學平衡方程表示為

(34)

由式(30)可以清楚地看出，式(34)中的每個方程都成立。節點各自由度對應的力和力矩狀態是平衡的，因此這種方法稱為力力矩平衡耦合方法。

對于線彈性靜力學問題，可以通過求解線性方程獲得模型的位移解，即

=

(35)

對于裂紋動態擴展問題，可采用中心差分法求解。系統待求解的運動方程為

(36)

4 算 例

4.1 帶初始裂紋板的裂紋擴展

本算例被廣泛用于評估裂紋擴展仿真算法的準確性。初始裂紋板模型及PD區域設置見圖6，板的材料參數為=72 GPa，=033，=2 440 kg/m，=135 J/m；作用在模型上下邊界的均布載荷為=12 MPa；將近場域大小設置為=3Δ。

圖 6 初始裂紋板模型,mm

計算得到不同時間步下的裂紋擴展路徑，見圖7。IGA-PD耦合模型模擬殼裂紋擴展是有效的。同時，根據裂紋擴展區域的分布，可選擇合適的區域作為PD區。對于裂紋擴展區域分布集中的問題，可以采用IGA-PD耦合算法提高裂紋擴展分析的計算效率。

圖 7 初始裂紋板的裂紋擴展路徑

4.2 平行預置裂紋板的裂紋擴展

為檢驗耦合算法在復雜裂紋條件下對裂紋擴展和合并的仿真能力，采用受位移載荷的預置平行裂紋板算例，見圖8。考慮裂紋角度=0和=45°這2種情況，幾何參數=10 mm，=10 mm。模型參數=30 GPa，=033，=2 700 kg/m，=30 J/m；將近場域大小設置為=3Δ。模型上下邊界給定0.04 mm位移，分2 000步線性加載，時間步長Δ=5×10s，則位移加載速度為0.4 m/s。PD區域根據預制裂紋分布進行布置。對于圖8中的傾斜裂紋，可設置包圍裂紋的區域為PD區域。裂紋擴展過程見圖9和10。2種角度的裂紋均在開始階段各自生長，內側相互靠近并逐漸合二為一，外側向兩側擴展至邊緣。裂紋擴展路徑與文獻[23]中的結果一致。

圖 8 平行預置裂紋板模型,mm

圖 9 α=0°平行裂紋板的裂紋擴展路徑

圖 10 α=45°平行裂紋板的裂紋擴展路徑

采用=0的模型研究不同PD區域大小對計算效率的影響。所用測試設備參數為CPU Intel i5-7400 3.00 GHz，RAM 32 GB 2 400 MHz，硬盤TOSHIBA DT01ACA100 7200 HDD 1 TB。以劃定的PD區域面積占整個模型總面積的比例為橫坐標，以采用耦合模型求解時間與采用純PD模型計算時間的比值為縱坐標，不同PD區域占比的計算時間與純PD區域的計算時間的相對用時對比見圖11。由此可知，隨著PD區域所占比例的增加，計算時間也增加。采用20%區域為PD區域，其余區域為IGA模型區域時，所需計算時間為純PD模型的25%，說明與PD方法相比，使用IGA-PD耦合方法能夠提高計算效率。

圖 11 不同PD區域時的計算時間對比

4.3 受內壓的圓柱殼

為檢驗耦合算法在一般殼結構斷裂仿真方面的性能，采用圓柱作為研究對象，模擬其在受內部壓強作用下的斷裂過程。帶預制裂紋的受內壓圓柱殼模型見圖12。材料參數=70 GPa，=03，=1.5 J/m；模型承受的內部壓強=1 MPa；設置近場域大小為=3Δ。

圖 12 帶預制裂紋的受內壓圓柱殼模型,mm

根據模型的對稱性，采用1/2幾何模型，軸向設置為67個控制點，周向設置35個控制點，得到2階張量積B樣條曲面。圓柱模型的裂紋擴展路徑見圖13。裂紋沿著與預置裂紋方向一致的方向擴展，并自始至終沿著與軸線平行的方向向兩側傳播。在一般殼結構中，所提出的耦合算法在精確幾何表示和斷裂仿真方面的性能均較好。

圖 13 圓柱模型的裂紋擴展路徑

5 結束語

針對殼結構中的裂紋擴展問題，建立等幾何分析(IGA)與近場動力學(PD)耦合的分析模型IGA-PD。在計算過程中，將模型域劃分為IGA域與PD域。IGA域采用基于Reissner-Mindlin(R-M)理論的退化殼模型，考慮橫向剪切的影響，保證幾何精確性，避免傳統有限元離散引入的幾何誤差，在法向向量計算等方面更加準確。PD域采用基于微梁鍵的非連續伽遼金(DG)有限元殼模型，實現對非連續場的計算分析。耦合過程基于力/力矩平衡原理，將等幾何剛度矩陣與近場動力學剛度矩陣耦合。數值試驗結果表明，該耦合模型可高效模擬板殼結構中的裂紋擴展，準確計算裂紋的開裂和交叉等現象。與PD模型相比，本文方法可有效提高計算效率，對于裂紋擴展區域相對集中的問題，可將計算時間縮短至原耗時的25%。