周倍合， 朱其志， 李惟簡

（1.河海大學巖石力學與堤壩工程教育部重點實驗室，南京 210098；2.河海大學江蘇省巖土工程技術工程研究中心，南京 210098）

固體材料破裂問題是固體力學研究的經(jīng)典難題之一［1］，根據(jù)研究手段不同現(xiàn)有研究方法大致分為兩類：物理試驗和數(shù)值模擬. 相比于試驗方法，數(shù)值模擬方法費用低、效率高，對計算模型的形狀和尺寸沒有特殊要求，邊界條件與荷載的設置更加自由，突破了試驗設備條件和試件制備方法的局限性. 目前，已出現(xiàn)多種數(shù)值方法，其中擴展有限元法（XFEM）和離散元法（DEM）是發(fā)展相對成熟并廣受關注的方法. XFEM方法在有限元方法（FEM）的基礎上引入額外的節(jié)點自由度和局部強化函數(shù)，可用于解決裂紋、孔洞、夾雜等間斷問題，無須重新劃分網(wǎng)格［3］，但不適用于模擬復雜交叉、分叉等多裂紋擴展貫通問題. 基于非連續(xù)性假設的DEM［6］方法可模擬裂紋的擴展貫通現(xiàn)象，但在反映材料連續(xù)區(qū)域的力學行為時存在明顯缺陷.

Silling 等［8］提出基于非局部理論的PD 理論，使用空間積分方程重構(gòu)固體力學運動方程，使方程在不連續(xù)處也有定義，實現(xiàn)在統(tǒng)一框架下自然地描述裂紋擴展過程. 早期的PD理論被稱為“鍵基”模型，存在泊松比固定的缺陷［10］. 為了克服泊松比限制問題并使PD理論與經(jīng)典連續(xù)介質(zhì)力學的聯(lián)系更加緊密，Silling等［9］提出了“態(tài)”的概念并發(fā)展了“態(tài)基”模型. 除態(tài)基PD模型外，學者們還提出多種改進的鍵基PD模型以解決泊松比限制問題，Han等［11］對這些模型進行了詳細介紹. 其中，Zhu等［12］提出了考慮鍵切向變形的雙參數(shù)鍵基PD模型（XPDM），在均勻變形條件下該模型可重現(xiàn)任意泊松比材料的變形，但在模擬非均勻變形問題存在不足. Diana等［13］提出的MPPD模型在XPDM模型的基礎上增加了局部轉(zhuǎn)動自由度，可考慮鍵的剛體轉(zhuǎn)動對自身切向變形的影響，并使用該方法成功模擬了非均勻變形問題.

使用顯式差分技術對運動方程進行時域積分是近場動力學的常規(guī)數(shù)值求解方法，對于靜態(tài)和準靜態(tài)問題，需要使用動態(tài)松弛技術對系統(tǒng)引入人工阻尼加快計算收斂速度［14］. 選擇合適的人工阻尼和時間步長對收斂速度和計算精度具有很大影響. 為解決這一問題，本文從能量角度出發(fā)，依據(jù)最小勢能原理推導了MPPD 模型物質(zhì)點系統(tǒng)的控制方程，構(gòu)造了系統(tǒng)的非局部剛度矩陣，實現(xiàn)用隱式方法求解靜態(tài)和準靜態(tài)問題. 使用該隱式方法對懸臂梁進行了彈性變形分析，并對雙懸臂梁試件中的Ⅰ型裂紋和L-形板中的Ⅰ-Ⅱ復合型裂紋擴展過程進行了模擬.

1 微極鍵基近場動力學模型

根據(jù)近場動力學方法的非局部特點，材料區(qū)域R 中的某個物質(zhì)點i 將與其近場域Hx（在二維和三維下分別為具有特征半徑δ 的圓和球體）中的其他質(zhì)點產(chǎn)生相互作用，這種相互作用通過“鍵”進行傳遞. 與只考慮軸向力的BB-PD模型不同，MPPD模型中考慮了鍵的切向變形和物質(zhì)點的旋轉(zhuǎn)，解除泊松比限制問題.可將MPPD模型中的鍵理想化為三個彈簧的組合，以鍵的初始狀態(tài)為參考，三個彈簧分為：約束物質(zhì)點沿著鍵的軸向和切向發(fā)生相對運動的彈簧以及約束物質(zhì)點發(fā)生相對旋轉(zhuǎn)的彈簧［13］，如圖1所示.

圖1 微極鍵基近場動力學模型示意圖Fig.1 Diagrammatic sketch of MPPD model

根據(jù)牛頓第二定律，t時刻物質(zhì)點i 的運動方程為：

其中：ρ 為物質(zhì)密度；J 為單位體積該物質(zhì)的轉(zhuǎn)動慣量；u?和θ?分別表示加速度和角加速度；f 和m 分別表示物質(zhì)點i 和j 之間的相互作用力和相互作用力矩；b 和c 分別表示物質(zhì)點i 處的體力密度和力矩密度；Vj為物質(zhì)點j 的體積. 當式（3）中兩個方程等號的左端項都為零時，物質(zhì)點i 處于平衡狀態(tài).

在小變形條件下，根據(jù)式（2）可得鍵的軸向伸長率l 可表示為：

式中：n∥=ξ ||ξ =(cos α,sin α)是初始狀態(tài)下鍵的單位方向向量；α 為初始狀下鍵的方位角. 規(guī)定逆時針旋轉(zhuǎn)方向為正，兩個質(zhì)點在t時刻的轉(zhuǎn)角分別為θi和θj，則小變形情況下鍵的切向變形可表示為：

式中：n⊥=(-sin α,cos α)，取n⊥的方向為切向變形的正方向. 以n∥和n⊥為正交基，可以建立鍵的局部坐標系統(tǒng)，如圖1（b）. 兩質(zhì)點的相對旋轉(zhuǎn)角為：

對于均質(zhì)各向同性線彈性材料，假設組成鍵的三個彈簧的勁度系數(shù)分別為kn、kt和k?，在t 時刻物質(zhì)點i 所受的彈簧力和力矩分別為：

式中：向量外積n∥×n⊥表示力矩的正方向. 根據(jù)上述條件，物質(zhì)點i 的應變能密度可表示為：

在均勻各向同性膨脹條件下γ=?=0，此時MPPD模型的應變能密度只與l 有關. 將此時MPPD模型的應變能密度與經(jīng)典連續(xù)介質(zhì)力學的應變能密度比較，可得到鍵參數(shù)kn與材料宏觀彈性參數(shù)之間的關系. 當kn已知，對于均勻的純剪切變形情況?=0，同樣可通過比較應變能密度可得到kt的表達式［13］.

平面應力情況下：

平面應變情況下：

在（9）和（10）式中：E 為楊氏模量；ν 為泊松比；h 為模型的厚度. 引入物質(zhì)點旋轉(zhuǎn)自由度可提高對非均勻變形問題模擬的準確性，平面問題中k?的表達式為［15-16］：

2 隱式求解方法

2.1 系統(tǒng)控制方程

數(shù)值求解時需要將模型離散化，若離散后任意質(zhì)點i 的近場域中有Mi個物質(zhì)點，式（8）的離散形式為：

式中：ξij= ||xj-xi為初始狀態(tài)時物質(zhì)點i 和j 之間的距離；?ij=θj-θi表示變形后物質(zhì)點i 和j 的相對轉(zhuǎn)角.υc(j)為物質(zhì)點j 的體積修正系數(shù)，表達式為［17］：

其中，r=Δx 2，Δx 為物質(zhì)點間距. 因為參數(shù)kn、kt和k?是基于具有完整近場范圍物質(zhì)點的應變能密度確定的，所以對于靠近自由表面或者材料的界面處的物質(zhì)點需要進行鍵參數(shù)修正. 本文使用能量法［17］計算表面效應修正因子Gij.

假設離散后模型中共有N 個物質(zhì)點，則物體的總勢能Etot(X)可表示為：

式中：X=(x′1,x′2,…,x′N)T表示物質(zhì)點系統(tǒng)的當前構(gòu)型；x′i(i=1,2,…,N)為變形后質(zhì)點i 的位置矢量；Vi表示質(zhì)點i 的體積；bi和ui分別表示質(zhì)點i 的體力密度矢量和位移矢量.

因為近場動力學中物質(zhì)點之間存在的非局部作用與原子尺度有限元方法［18］中的多體相互作用類似，根據(jù)最小勢能原理，當離散后的物質(zhì)點系統(tǒng)總勢能最小時物體處于平衡狀態(tài). 假定X=X?時物質(zhì)點系統(tǒng)的總勢能最小，則

將式（16）代入式（15）得到整個系統(tǒng)的控制方程：

其中，位移u=X-X0，K 為系統(tǒng)的非局部剛度矩陣：

P 為系統(tǒng)的不平衡力，可表示為：

對于非線性系統(tǒng)，需要對式（17）進行迭代計算，每個迭代步中更新剛度矩陣和不平衡力，直到系統(tǒng)不平衡力P 等于零.

2.2 單元剛度矩陣和不平衡力

將物質(zhì)點看作有限元節(jié)點時，式（17）與經(jīng)典有限元方法的控制方程一致. 有限元方法常用的單元形式只能描述鄰近節(jié)點之間的局部作用，不適用于考慮非局部作用PD模型. 為了描述這種非局部作用，可將中心質(zhì)點i 與其近場域中所有質(zhì)點組成的集合作為一個單元，如圖2所示. 單元中，中心點與近場域中所有非中心點間存在相互作用，但非中心質(zhì)點之間互不影響. 因此，在單元剛度矩陣和不平衡力中，只有與中心點有關聯(lián)的位置上元素不為零，其余的元素均為零. 物質(zhì)點i 處對應單元的剛度矩陣和不平衡力分別為：

圖2 MPPD-FEM模型中的一個單元Fig.2 Illustration of the element used in MPPD-FEM model

其中，j=1,2,…,Mi，Mi表示物質(zhì)點i 近場域中包含其他物質(zhì)點的個數(shù). 在MPPD模型中每個物質(zhì)點具有三個自由度，物質(zhì)點i 對應單元的剛度矩陣為3(Mi+1)階的方陣.

2.3 損傷的描述

近場動力學中，材料損傷可通過消除物質(zhì)點之間的相互作用來描述. 在微觀彈脆性（PMB）模型中，物質(zhì)點之間鍵的伸長率超過臨界值lc時，鍵失效［19］. 引入標量特征函數(shù)λ(η,ξ)來描述鍵的失效情況：

由于MPPD模型解除了泊松比限制，不同泊松比時的臨界伸長率可按以下方法計算［20］：

其中，Gc為斷裂力學中的能量釋放率；β=3δ 4π；β′=0.238 73δ；μ 和κ 分別為剪切模量和體積模量：

物質(zhì)點的局部損傷定義為該物質(zhì)點近場域中失效鍵的數(shù)量與初始狀態(tài)時鍵總數(shù)的加權(quán)比值：

使用微觀彈脆性（PMB）模型能夠有效模擬巖石、混凝土等脆性材料的破壞問題. 為簡化問題，本文基于鍵基模型的微觀脆性材料進行研究，使用最大伸長率準則判斷MPPD模型中鍵的軸向彈簧kn是否發(fā)生破壞，并假設軸向彈簧破壞時切向彈簧和旋轉(zhuǎn)彈簧也一同破壞，對應物質(zhì)點之間的相互作用失效.

2.4 隱式方法的實現(xiàn)流程

微觀彈脆性材料的試件在出現(xiàn)損傷前為線性系統(tǒng)，荷載作用下發(fā)生線彈性變形. 隨著荷載增加，物質(zhì)點系統(tǒng)中有達到臨界伸長率的鍵失效時，材料開始出現(xiàn)損傷. 繼續(xù)加載過程，失效鍵的數(shù)量不斷增加，材料中損傷累積并導致力學性能逐漸衰減. 因此，當計算過程中有新的鍵失效時，需要及時更新整體剛度矩陣，將失效鍵對應的貢獻值從中減去. 為了準確模擬材料的損傷變化過程并且兼顧計算效率，參考文獻［21］提出的C方案，每個迭代步中只斷開達到臨界伸長率且伸長率最大的Nb個鍵，然后更新剛度矩陣. 當某個迭代步中沒有新的鍵失效時，增加荷載然后再進行修正剛度矩陣的迭代過程，直到荷載大小達到停止條件或試件完全破壞時結(jié)束程序. 隱式算法的程序?qū)崿F(xiàn)過程如圖3所示，本文在計算過程中取Nb=4 .

圖3 隱式求解方法的流程圖Fig.3 Flow chart of the implicit solution method

3 算例與分析

3.1 彈性變形分析

如圖4 所示，平面應力狀態(tài)的懸臂梁，跨長為L=1 m，高度H=0.3 m，彈性模量E=200 GPa，泊松比ν=0.2. 梁的左端完全固定，右端自由并受到豎直向下的均勻分布的荷載，荷載合力F=1 kN. 用均勻分布的物質(zhì)點對模型進行離散，物質(zhì)點間距Δx=0.01 m，近場域半徑δ=3.015 Δx .

圖5（a）為使用MPPD方法模擬得到的沿x 和y 方向位移分量分布云圖，F(xiàn)EM 法所得結(jié)果如圖5（b）所示.對比可知，兩種不同方法的模擬結(jié)果十分接近，位移場分布規(guī)律相同，但在數(shù)值上存在細微差別. MPPD方法計算的沿x 和y 方向的位移分量最大值都略小于有限元的結(jié)果，但相對于有限元的誤差不超過1.0% .

圖4 懸臂梁示意圖Fig.4 Geometry and boundary conditions of a cantilever beam

3.2 裂紋擴展模擬

3.2.1 雙懸臂梁中Ⅰ型裂紋擴展模擬 平面應力狀態(tài)雙懸臂梁的幾何信息與邊界條件如圖6所示，厚度取h=10 mm，彈性模量E=200 GPa，泊松比ν=0.2，能量釋放率Gc=962.8 J/m2. 預留初始裂紋位于試件的水平中心線上，長度為100 mm. 用均勻分布的物質(zhì)點對模型進行離散，物質(zhì)點間距Δx=0.5 mm，近場域半徑δ=3.015 Δx. 雙懸臂梁試件中裂紋先沿著水平方向進行開展，當裂紋接近試件固定端時將會發(fā)生彎折，如圖7所示. 裂紋沿直線擴展階段的位移-荷載曲線見圖8，其結(jié)果與Abaqus中基于虛擬裂紋閉合技術的有限元分析結(jié)果相吻合.

圖5 懸臂梁模型的位移分量分布云圖（單位：m）Fig.5 Node displacements obtained by MPPD and FEM

圖6 雙懸臂梁幾何示意圖Fig.6 Geometry and boundary conditions of a double cantilever beam

圖7 雙懸臂梁試件中的裂紋路徑Fig.7 Crack paths of double cantilever beam

圖8 雙懸臂梁裂紋開展的位移-荷載曲線Fig.8 Displacement-load curves of double cantilever beam crack propagation

3.2.2 L-形板中裂紋擴展模擬 普通混凝土材料的L-形板試件幾何信息與邊界條件如圖9（a）所示，厚度為h=100 mm，彈性模量為E=25.85 GPa，泊松比ν=0.18，能量釋放率Gc=65 J/m2，按平面應力問題處理.模型用均勻分布的物質(zhì)點離散，物質(zhì)點間距Δx=2.5 mm. 為了準確模擬L-形板中混合型裂紋的開展情況，適當增大了近場域的半徑，取δ=5.015 Δx. 圖10展示了不同位移條件時的模擬結(jié)果，由圖可知，裂紋從L-形板的角點處出現(xiàn)，初始階段與水平線成一定角度（約為26°），隨著加載位移增大，裂紋開展方向逐漸趨于水平并向試件的右邊界擴展. 模擬結(jié)果與圖9（b）中試驗結(jié)果［22］對比具有較好的一致性.

圖9 L-形板示意圖和試驗觀測的裂紋形狀Fig.9 Sketch map of L-shape structure test and the experimentally observed crack path

圖10 L-形板中裂紋擴展模擬結(jié)果Fig.10 Crack propagation mode in L-shape structure test

4 結(jié)論

1）本文從能量角度出發(fā)，依據(jù)最小勢能原理重新推導了MPPD模型物質(zhì)點系統(tǒng)的控制方程，提出一種適用于MPPD模型的隱式求解方法，用于提高靜態(tài)和準靜態(tài)問題的求解效率. 提出一種可以描述物質(zhì)點之間非局部作用的單元形式，不再以單獨的近場動力學鍵為單元，可提高組裝整體剛度矩陣的效率.

2）分別使用隱式求解格式的MPPD方法和有限元方法對非均勻變形的懸臂梁進行彈性響應分析，得到物質(zhì)點沿x 和y 方向位移分量的分布云圖. 不同方法模擬得到的位移場分布規(guī)律相同，在數(shù)值上存在微小差別，但滿足計算誤差要求，驗證了MPPD模型隱式求解方法的有效性.

3）使用按隱式方法求解的MPPD模型對雙懸臂梁試件中的裂紋擴展現(xiàn)象進行模擬，裂紋沿預期方向穩(wěn)定開展，并且裂紋水平直線階段的位移-荷載曲線與Abaqus中基于虛擬裂紋閉合技術的有限元分析結(jié)果吻合良好. 此外，還對普通混凝土材料的L-形板進行了模擬，所得混合型裂紋的開展路徑與試驗結(jié)果十分接近. 通過以上兩種算例驗證了MPPD 模型隱式求解方法用于模擬靜態(tài)和準靜態(tài)條件下裂紋擴展問題的可靠性.