近場動力學與有限元混合模型研究

錢 偉，范存新，沈 峰，夏益兵

(蘇州科技大學 土木工程學院，江蘇 蘇州 215000)

以有限單元法為代表的基于偏微分方程的相關數值方法已經得到了廣泛的運用，但在模擬裂紋擴展等不連續問題時，基于連續性假設的有限單元法無法自發模擬裂紋擴展，國內外學者通過設置界面單元或采用網格重劃分技術來處理，但仍存在網格依賴性的問題[1]。

2000年，美國Sandia國家實驗室的Silling博士提出了基于非局部作用思想的近場動力學方法[2-3]，在國際上引起了廣泛關注，從根本上解決了連續介質力學在模擬裂紋路徑等不連續問題時偏導數不存在的問題，逐漸成為計算力學與工程仿真及相關領域研究熱點[4]，在斷裂破壞問題中得到了廣泛應用[5-8]。然而，近場動力學計算效率相較有限元而言過低。為此，國內外學者嘗試將近場動力學與有限元結合起來，充分發揮兩者優勢。Macek和Silling[9]將研究對象劃分為PD子域與FE子域以及重疊域，重疊域采用有限元實體模型，將重疊域中的PD模型以及PD子域視作桿單元，通過鑲嵌單元技術，實現了對PD與FEM的混合建模；Liu和Hong[10]在近場動力學子域與有限元子域設置界面單元，在界面單元內加入一定的物質點，計算耦合力并將耦合力通過兩種方式分配到界面單元的結點上，以此來實現兩種方法的混合建模；Seleson等[11-12]通過在[0，1]區間內變化的混合函數將局部作用區域的應變能密度與PD應變能密度混合實現了兩區域的平滑過渡；Ren等[13]利用不連續的Galerkin法建立了經典近場動力學控制方程，并給出了相關算例。上述混合模型均在顯式體系下求解，處理靜力問題時需引入人工阻尼，阻尼的大小影響迭代收斂速度，進而導致計算效率不足。此外，上述文章均局限于通過低階單元與近場動力學方法實現混合建模，缺乏對高階單元與近場動力學混合建模研究。

在此基礎上，本文所建立的模型將研究對象劃分為PD子域與FE子域，在裂紋出現的區域，采用近場動力學模型，其他區域采用等參元模型。通過桿單元連接PD物質點與等參元結點來實現混合建模。該模型無需引入人工阻尼，提高了計算效率。同時，該混合模型在FE子域采用八結點等參元，提高了計算精度，最后，通過對懸臂梁彈性變形及含裂紋正方形板破壞過程的模擬，證明了該混合模型的有效性。

1 近場動力學簡介

1.1 近場動力學方法

如圖1所示，設在任意時刻t,任意空間R內任一物質點x與其鄰域一定范圍δ內的其他物質點(x′∈R:‖x′-x‖≤δ)存在單位相互作用力f，根據牛頓第二定律，可得[14]：

圖1 物質點間的相互作用Fig. 1 Interaction of material points

(1)

令式(1)中x′-x=ξ，u′-u=η，ξ為物質點間相對位置，η為物質點間相對位移。

式中：Hx為物質點x的鄰域范圍，δ為鄰域范圍尺寸，ρ為物質點材料密度，u為物質點的位移，b為外力密度。

基于保守場的定義和性質，一定存在一個標量函數w(η,ξ)(物質點對勢能密度)，使得：

(2)

物質點間的作用類似一個中心彈簧，則

(3)

(4)

μ為一標量函數，用來表征鍵的斷裂：

(5)

式中:s0為臨界伸長率。當伸長率小于臨界伸長率s0時，μ(t,ξ)=1，鍵未發生斷裂，否則，μ(t,ξ)=0，鍵斷裂。這里考慮材料的拉壓異性，可以通過材料的抗拉強度、抗壓強度和彈性模量E去表征：

(6)

在近場動力學理論中，統一定義局部損傷

(7)

式中：0≤φ≤1，φ=0表示材料未損傷，φ=1表示該點完全損傷。

1.2 改進的PMB模型

式(3)中的微觀模量函數c(ξ)可表示為：

c(ξ)=c(0,δ)g(ξ,δ)，

(8)

式中：c(0,δ)為集中函數，g(ξ,δ)為核函數，表示遠程力大小隨物質點間距變化的規律。在改進的PMB模型中，考慮遠程力對微觀模量的影響，取核函數[16-17]為：

(9)

根據近場動力學應變能密度與連續介質力學應變能密度相等的原則，可以得到平面應力狀態下的集中函數為：

(10)

2 有限元簡介

2.1 有限元方程

有限單元法靜力問題求解方程為

KU=R，

(11)

式中：K為整體剛度矩陣，U為整體位移列陣，R為整體等效結點荷載列陣。K和R可由式(12)得出

(12)

式中：Ce為選擇矩陣，k為單元剛度矩陣，Re為等效荷載列陣，B為應變轉換矩陣，D為彈性矩陣。

2.2 八結點等參元

相較于矩形單元，等參元精度高，且適用于復雜的曲線邊界與曲面邊界，因此得到了廣泛的應用。這里在有限元子域采用八結點等參元。

圖2為四邊形單元在整體坐標系xy下的單元結點分布，在每個四邊形單元上建立局部坐標系ξη，如圖3所示，通過坐標變換，將每個四邊形單元映射到標準正方形單元，建立了兩種單元的對應關系。

圖2 四邊形單元Fig. 2 Quadrilateral element

圖3 正方形單元Fig. 3 Square element

設單元中任意一點的位移是u，v，單元結點位移為ui，vi(i=1，8)，其位移模式為

(13)

坐標變換為：

(14),

式中，Ni=(i=1,8)為八結點等參元形函數，其表達式為：

(15)

式中，ξ,η是定義在標準單元上的局部坐標，ξi，ηi(i=1,2,3,4)分別代表標準單元4個結點的局部坐標值。

3 近場動力學與有限元耦合方案

如圖4所示，將幾何模型劃分為PD子域與FE子域，在PD子域采用近場動力學建模，在FE子域采用八結點等參元建模。在兩類子域的交界面上，如圖5所示，交界面上的每個有限元結點通過桿單元與其近場范圍內的物質點相連接[18-19]，其中，PD子域物質點對間相互作用可視為桿單元。為求解近場動力學方程，將研究對象離散為一系列帶有物性信息的物質點，近場動力學積分方程轉化為對有限個物質點的求和，即：

圖4 有限元(FE)與近場動力學(PD)子域Fig. 4 FE and PD region

圖5 FE結點與PD物質點連接示意圖Fig. 5 Interactions between FE nodes and PD nodes through trusses

(16)

對于靜力問題，令加速度為0，可得近場動力學的平衡方程：

(17)

根據公式(1)，可將PD方程改寫成矩陣形式

(18)

(19)

桿單元的剛度貢獻矩陣為：

(20)

通過求解式(12)、式(19)、式(20)，分別得到了等參元、物質點對間、桿單元的剛度貢獻矩陣，其中，物質點對間相互作用視為桿單元，通過對剛度貢獻矩陣的集成，形成整體剛度矩陣。最后根據有限元靜力方程求解位移，實現了近場動力學與有限元的混合建模。

4 數值算例

4.1 懸臂梁的彈性變形

懸臂梁幾何尺寸及子域劃分如圖6所示[19]，跨長1 000 mm，截面高為200 mm，彈性模量E為100 GPa，泊松比ν=1/3，物質點間距取為2.5 mm，近場范圍尺寸取δ=4Δx，有限元網格為10 mm×10 mm的八結點等參單元，右端受大小為1 050 kN/m的均布荷載，探討不同數值方法對精度的影響。

圖6 懸臂梁幾何模型Fig. 6 Geometric model of cantilever beam

表1給出了采用不同數值方法得到的最大水平位移計算結果。由表1可得，采用近場動力學方法計算的最大水平位移與解析解的相對誤差為2.19%，采用四結點混合模型時相對誤差為0.67%，采用八結點混合模型的誤差為0.28%，混合模型的精度比近場動力學的計算精度高，且本文中建立的八結點混合模型高于四結點混合模型的計算精度。證明了本文中提出的混合模型的精確性。其中，當不設置FE子域時，可根據有限元靜力求解方程，得到近場動力學計算結果。

表1 最大水平位移計算結果

4.2 含I型裂紋板受拉破壞分析

考慮如圖7所示含I型裂紋正方形板，模型尺寸及材料參數如下，邊長為50 mm，在板中間預制一條長為10 mm的裂紋，彈性模量為30 GPa，ft為2.01 MPa，fc為20.1 MPa，泊松比為1/3。將正方形板劃分為兩個有限元子域與一個近場動力學子域，物質點間距取為0.5 mm，近場范圍尺寸取δ=4Δx，有限元網格為2 mm×2 mm八結點等參元，采用位移加載，每一步位移增量為1.0×10-8m。

圖7 正方形板幾何模型Fig. 7 Geometric model of square plate

混合模型計算得到的正方形板裂紋擴展過程如圖8所示，當加載到57步時(此時位移荷載為5.7×10-7m),預制裂紋的裂尖出現損傷；隨著荷載的進一步增大，裂紋發生擴展，如圖8(b)(c)所示；當加載至92步時，(此時位移荷載為9.2×10-7m)，裂紋貫穿整個正方形板，構件發生破壞。

圖8 裂紋擴展示意圖Fig. 8 Crack propagation process

4.3 多裂紋板受拉破壞分析

考慮如圖9所示的含兩條斜裂紋的正方形板，模型尺寸及材料參數如下，板的邊長為50 mm，初始裂紋長度為10 mm，b為10 mm，彈性模量為30 GPa，ft為2.01 MPa，fc為20.1 MPa，泊松比為1/3。將正方形板劃分為兩個有限元子域與一個近場動力學子域，物質點間距取為0.5 mm,近場范圍尺寸取δ=4Δx，有限元網格為2 mm×2 mm的八結點等參單元，采用位移加載，每一步位移增量為1.5×10-8m。

圖9 正方形板幾何模型Fig. 9 Geometric model of square plate

混合模型計算得到的正方形板裂紋擴展過程如圖10所示，加載到39步時(此時位移荷載為5.85×10-7m)，裂紋裂尖出現損傷；隨著荷載的進一步增大，裂紋發生擴展并交匯，如圖10(b)(c)所示；當加載至63步時，(此時位移荷載為9.45×10-7m)，裂紋貫穿整個正方形板，構件發生破壞。圖10(e)為多維虛內鍵計算結果[20]。對比可得，混合模型與虛內鍵計算結果基本吻合。

圖10 裂紋擴展示意圖Fig. 10 Crack propagation process

5 結 論

PD理論通過求解積分方程模擬材料斷裂破壞行為，在分析裂紋擴展等不連續問題時具有顯著優勢，然而，PD理論計算效率相較FEM而言過低。為兼顧兩者優勢，采用近場動力學與有限元混合建模的方法，建立了新的混合模型，該模型將研究對象劃分為PD子域與FE子域，在PD子域采用近場動力學建模型，FE子域采用等參元模型，通過桿單元連接PD物質點與等參單元結點來實現混合建模。此混合模型在求解處理靜力問題時無需引入人工阻尼，采取類似有限元的方法，對剛度集成，形成整體剛度矩陣，然后根據有限元支配方程求解靜力問題，提高了計算效率。同時在FE子域采用八結點等參元建模，提高了混合模型的計算精度。最后，采用所建立的混合模型計算分析了懸臂梁的彈性變形和含裂紋正方形板的破壞過程，取得了較好的結果，為斷裂破壞問題的解決提供了一種新思路。