李樹忱，馬鵬飛，王修偉，劉祥坤

(山東大學巖土與結(jié)構(gòu)工程研究中心，山東，濟南 250061)

固體材料的力學行為一直是工程結(jié)構(gòu)領域關注的熱點，探究其失穩(wěn)破壞機理對工程防災減災具有重要意義。多年來，國內(nèi)外學者針對此問題已取得了豐碩研究成果[1-4]，采用理論方法探究材料破壞機制時通常需要一定簡化，無法很好地反映真實的情況，利用試驗手段存在成本較高、時間周期長等困難。隨著計算機技術發(fā)展，數(shù)值模擬手段越來越多的被用來研究材料的力學特性。

有限元法[5-7]、離散元法[8-9]、無網(wǎng)格方法[10]等數(shù)值方法的相繼出現(xiàn)為探究材料力學行為提供了成本低且較快捷的方式，在模擬裂紋擴展時微分控制方程會遇到尖端奇異性等難題。擴展有限元理論運用預先設定裂紋擴展方向及長度來模擬材料破壞現(xiàn)象[11-12]。分子動力學方法可較好的體現(xiàn)結(jié)構(gòu)微觀信息，但是在模擬裂紋擴展等尺度時會遇到計算效率低、耗費時間長等不足[13]。

近場動力學方法是采用非局部作用特性來描述材料力學行為的理論[14]。該方法核心為求解空間積分方程，因此可較好的避免尖端奇異性等難題。近場動力學理論最早由SILLING 等[15]提出，黃丹等[16-17]首次將該理論引入國內(nèi)，并針對混凝土等固體材料變形破壞開展研究。ZHOU 等[18]和谷新保等[19]利用近場動力學方法模擬了巖石類材料的裂紋擴展過程。WANG 等[20-21]建立了共軛鍵基近場動力學模型以突破泊松比的限制。ZHU等[22]提出了具有鍵旋轉(zhuǎn)效應的近場動力學公式并對脆性材料變形行了模擬。王超等[23]在模擬潛艇破冰上浮的現(xiàn)象中運用了近場動力學方法。黃小華等[24]提出了雙參數(shù)模型并對沖擊荷載對破壞的影響進行了分析。HAN 等[25]利用連續(xù)介質(zhì)力學理論與近場動力學耦合的方式對材料變形進行模擬。

近場動力學中以鍵為基礎的微彈性模型具有廣泛的應用[26]，微彈性模型控制方程為積分方程，在計算時不涉及應力及應變，但經(jīng)典固體材料力學特性分析時應力與應變更為直觀，并且大多破壞準則為基于應力的準則，在近場動力學理論中無法直接應用，文獻[27]中提出利用力通量的幾何等價關系計算應力，但其計算過程復雜且計算量相對較大。

本文結(jié)合前人研究基礎，在近場動力學中引入非局部微分算子求解理論，建立微彈性應力分析模型。為驗證有效性，利用擴展后的方法對幾種二維平面問題進行模擬，同時將計算值與經(jīng)典理論解對比，并且對粒子離散間距、泰勒項數(shù)及權函數(shù)的數(shù)值收斂性進行研究。結(jié)果表明：本文提出的方法可較好的計算結(jié)構(gòu)處的應力值并進行相應的分析。

1 基本理論

1.1 近場動力學理論

近場動力學理論與傳統(tǒng)方法相比最大的不同為非局部作用。如圖1 所示非局部作用是指在一定范圍內(nèi)存在相互作用力，x′為x作 用域 δ內(nèi)具有物理意義的物質(zhì)點，在初始參考坐標系中兩點的相對位置為x′-x=ξ ，并且x在t時刻的控制方程為：

式中：u(x,t)與u(x′,t) 分 別為x與x′變形時的位移，u(x′,t)-u(x,t) =η為變形后的相對位移，因此，此時x與x′相對位置可表示 ξ+η ； ρ(x) 與u¨(x,t)分別為密度與加速度；f為x與x′之間相互作用力；b(x,t)為t時刻體積力密度。

1.2 微彈性模型

在微彈性模型中相互作用力為相對位置 ξ與位移 η的函數(shù)，同時x與x′之間f大小相等方向相反，因此，式(1)可重新表示為：

式中，1/2 為每個相互作用鍵中的能量一半屬于物質(zhì)點x，其中的微勢能ω可寫為：

2 非局部微分算子

2.1 積分方程

應力-應變在經(jīng)典Navier 平衡方程中為位移的微分函數(shù)，而非局部微分算子試圖尋找其積分表達以適用于近場動力學模型，本節(jié)以二維為例推導微分算子的形式。如圖1 所示，x與x′相互用點并且相對位置x′-x=ξ ，可考慮在x處進行Taylor級數(shù)展開：

2.2 應力-應變表達

由近場動力學與經(jīng)典Navier 平衡方程對應關系[28]，這里將近場動力學中某點應力分量表示為：

本節(jié)對完整平板、含圓孔平板及含裂隙平板在外荷載作用下的變形進行模擬，并將應力分布結(jié)果與經(jīng)典理論解對比以驗證本文方法的有效性?，F(xiàn)在考慮一個各向同性的平板在兩端受到均勻荷載其尺寸及荷載分布如圖2 所示。

平板尺寸1 m×0.5 m，采用平面應力條件計算，彈性模量E=18 GPa ， 泊松比ν=1/3 ，密度ρ=2730 kg/m3，兩側(cè)荷載為p0=15 MPa，以中心為坐標原點應力分布的理論解為：

3 數(shù)值驗證

3.1 平板拉伸

計算時將模型離散為 200×100節(jié)點，間距取Δx=0.005 m，相互作用半徑 δ=3.015Δx同時采用穩(wěn)態(tài)求解方案，計算結(jié)果對比如圖3 所示。

由圖3 可知利用本文方法所計算的應力結(jié)果與理論解相比在模型中間區(qū)域吻合良好，但在邊界及角落處存在較大誤差，主要原因是邊界附近的節(jié)點鄰域不完整導致模量值失真從而求解位移場出現(xiàn)偏差，因此，基于位移求解的應力場在邊界區(qū)域同樣會出現(xiàn)誤差。

3.2 含圓孔平板

為了進一步驗證本文方法的有效性，3.2 節(jié)對具有高次精確解的含圓孔平板在荷載下變形進行模擬并將計算結(jié)果與經(jīng)典理論解對比?，F(xiàn)在考慮一個中間含有孔的方板在兩端受到均勻荷載，尺寸及荷載分布如圖4 所示。

平板尺寸0.5 m×0.5 m 并且中間含R=0.025 m的圓孔，采用平面應力計算，彈性模量E=15 GPa，泊松比ν=1/3 ，密度 ρ=2650 kg/m3，兩側(cè)荷載為p0=10 MPa，以中心為極坐標原點應力分布的經(jīng)典理論解為：

由圖5 可知平板的孔口區(qū)域表現(xiàn)出明顯的應力集中現(xiàn)象，孔口應力遠大于距孔口較遠處的應力。同時最大和最小的應力均出現(xiàn)在孔邊上，由開孔所引發(fā)的引力擾動也主要發(fā)生在1.5 倍圓孔直徑的范圍之內(nèi)，從定性的方面驗證本文方法的有效性。由直角坐標系與極坐標系的對應關系，可分別選取豎直方向、水平方向與環(huán)繞孔口區(qū)域的計算值與理論值比較，從定量的角度驗證有效性，對比結(jié)果如圖6 所示。

如圖6 所示 σxx和 σyy沿豎直與水平方向?qū)Ρ龋?σxy沿半徑 1.4R的環(huán)向進行對比。利用本文方法所得到的結(jié)果與理論解的整體偏差很小，誤差僅出現(xiàn)在靠近圓孔邊緣附近。理論解圓孔附近應力集中系數(shù)為3，本文方法所預測的應力集中系數(shù)為2.7 左右，逐漸遠離孔口區(qū)域的計算結(jié)果收斂于理論解，靠近模型邊界時會發(fā)生略微的失真但整體吻合較好，因此，可以驗證本文方法的有效性。

3.3 含裂隙平板

實際結(jié)構(gòu)破壞時往往會伴隨裂紋的萌生與擴展，探究裂紋周圍應力場分布對裂紋的發(fā)展具有指導意義，因此，本節(jié)利用該方法計算裂紋尖端應力場并將結(jié)果與理論解對比驗證方法的有效性。

現(xiàn)在考慮一個中間含有裂紋的方板在四周受到均勻荷載，尺寸及荷載分布如圖7 所示。平板尺寸0.5 m×0.5 m 并且中間含L=0.1 m的裂紋，采用平面應力進行計算，彈性模量E=15 GPa，泊松比v=1/3 ， 密度 ρ=2600 kg/m3，四周施加的荷載為p0=1 MPa，該問題的經(jīng)典斷裂力學解為：

式中：a為預制裂隙一半的長度；KI為應力強度因子，計算得KI=0.396 MPa·m1/2； θ與r為裂紋尖端附近點的角度與半徑。計算時將模型離散為150×150 節(jié) 點，相互作用半徑 δ=3.015Δx，其余計算參數(shù)與3.2 節(jié)相同應力計算結(jié)果如圖8 所示。

由圖8 可知在應力作用下平板裂隙處存在明顯的應力集中現(xiàn)象，而且隨著與裂隙距離的增加逐漸收斂于邊界值，從定性方面驗證方法有效性。

由圖9 為理論解的局部對比，本文方法計算的結(jié)果與理論解相比在裂紋尖端處應力分布與最大值基本一致，而遠離尖端處時計算解則會與理論解產(chǎn)生偏差并逐漸收斂于邊界條件值。偏差產(chǎn)生的原因是經(jīng)典斷裂力學解的適用范圍為裂紋尖端處，并且有模型無限大的前提假定，模型尺寸相比裂紋不夠無限大導致理論解與計算結(jié)果不同，并且由于在力學分析中往往關注裂紋尖端處的應力狀況，因此，可以驗證本文方法的有效性。

3.4 裂紋擴展過程

前面幾節(jié)所計算應力均為穩(wěn)定狀態(tài)下的結(jié)果，為了驗證本文方法同樣適用于動態(tài)裂紋擴展，本節(jié)將探究含預制裂紋材料在壓縮破壞下的應力變化過程。計算模型在圖7 的基礎上進行修改，考慮一個中間含有裂紋板在豎向受持續(xù)位移邊界條件，平板尺寸0.5 m×1.0 m，裂隙長度 0.28 m，裂隙與豎直方向的為夾角 30°，彈性模量取E=4.71 GPa，密度取 ρ=2600 kg/m3，臨界伸長率取s0=0.0318，在持續(xù)位移加載下發(fā)生逐漸破壞，損傷及最大主應力分布過程如圖10 所示。

由圖10 可知，預制裂隙在豎向持續(xù)加載的條件下發(fā)生漸進破壞，裂紋擴展的方向與加載方向相同，其中最大主應力出現(xiàn)正負交替分布的情況(拉為正壓為負)。在損傷初期階段拉應力出現(xiàn)在裂紋端部靠近加載端的部位，隨著荷載的不斷增加拉應力不斷變大，裂紋出現(xiàn)萌生擴展，拉應力持續(xù)出現(xiàn)在裂紋尖端處促使裂紋向最大主應力方向擴展，這過程與經(jīng)典分析一致，可驗證方法有效性。

4 數(shù)值收斂性分析

4.1 離散間距影響

本節(jié)開展非局部微分算子離散間距、泰勒項數(shù)及權函數(shù)收斂性研究，探究不同因素對結(jié)果的影響。以3.2 節(jié)含圓孔平板為例，將模型分別離散為 50×50 、 80×80 及 100×100節(jié)點，其余計算參數(shù)與3.2 節(jié)相同，x=0.0方向結(jié)果對比如圖11 所示。

由圖11 可知，隨著離散節(jié)點的增加粒子離散間距隨之減小，節(jié)點更加密集，在x=0.0方向上σxx與 σyy的最終計算結(jié)果更趨近于理論解，這一結(jié)論同樣適用于y=0.0方向，在此不再贅述。由于增加節(jié)點數(shù)量會成倍的提高計算效率，因此，在實際計算時應該合理選擇離散間距兼顧效率與準確度。

4.2 泰勒項數(shù)影響

本節(jié)計算條件與3.2 節(jié)一致，其中泰勒項數(shù)N分別取2、3 及5，x=0.0方向?qū)Ρ热鐖D12 所示。

由圖12 可知，隨著泰勒項數(shù)的增加，對x=0.0 方 向上 σxx與 σyy最終計算結(jié)果的影響并不顯著。泰勒項數(shù)的增加會導致計算效率降低，因此，在實際計算時應該同樣應該合理選擇確保效率。

4.3 權函數(shù)影響

本節(jié)計算條件與3.2 節(jié)一致，其中權函數(shù)除了指數(shù)分布外，同時引入如下所示的Gauss 函數(shù)、三次樣條及四次樣條權函數(shù)，其中d=|ξ|/δ，探究權函數(shù)對結(jié)果的影響，結(jié)果對比如圖13 所示。

由圖13 可知，Gauss 權函數(shù)計算結(jié)果與理論解偏差最小，該結(jié)論在 σyy邊界區(qū)域較為明顯，計算效果比樣條權函數(shù)好，因此，在計算時可考慮采用Gauss 權函數(shù)。

5 結(jié)論

本文建立了基于非局部微分算子的近場動力學應力分析模型并對其有效性及收斂性進行了研究，取得以下結(jié)論：

(1) 基于非局部微分算子的近場動力學應力分析模型可以較好計算固體材料的應力分布狀況，適用于完整、非完整材料，為經(jīng)典近場動力學方法提供了應力分析思路。

(2) 離散間距及權函數(shù)對數(shù)值收斂結(jié)果具有顯著影響，其中Gauss 權函數(shù)具有較好地收斂效果，同時計算精度隨著離散間距的減小而提高但計算效率會降低，實際計算時應合理選擇離散間距。