馬鵬飛，李樹忱，王修偉，王曼靈

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

固體材料的力學特性始終是工程結構領域關注的熱點，探究其失穩(wěn)機制對工程防災減災領域意義重大[1-6]，采用理論方法探究裂紋擴展機制時通常需要大量不能真實反映工程問題，而利用實驗方式時存在成本高、時間長等難題。隨著計算機技術的發(fā)展，數(shù)值模擬方法越來越多的被用來探究結構變形及破壞問題。

有限元法[7-9]、離散元法[10-11]、無網(wǎng)格方法[12-13]等數(shù)值方法的出現(xiàn)為探究結構變形提供了低廉、快捷的途徑，但在模擬材料裂紋擴展時會遇到尖端奇異性及裂紋成核等難題。擴展有限元理論采用預先設定裂紋的擴展方向及長度來解決結構裂紋擴展問題[14-15]。分子動力學方法可以較好的反應結構微觀尺度，但是在處理宏觀尺度問題時遇到計算效率低、計算時間長等難題[16]。

近場動力學是利用非局部作用思想來描述物質力學行為的方法[17]。該理論求解空間積分型運動方程而非微分方程可以很好的避免裂紋尖端奇異性等問題。近場動力學最早由Silling 等[18]提出，黃丹等[19-20]將其引入國內，針對混凝土變形破壞、巴西劈裂等問題進行了模擬。Zhou 等[21]利用近場動力學理論模擬了巖石類材料的裂紋擴展問題。王超等[22]在冰槳接觸及模擬潛艇破冰上浮問題中應用了近場動力學。Zhu 等[23]運用該方法對不含不同預制傾角的巖石試件破壞過程進行了模擬。王超聰?shù)萚24]采用近場動力學理論模擬了熱防護材料燒蝕過程。黃小華等[25]提出了單鍵雙參數(shù)近場動力學模型并對沖擊荷載對破壞的影響進行了分析。

以鍵為基礎的微彈性模型(PMB)思路清晰應用廣泛[26]。由于該模型過度簡化，使得利用該模型時泊松比為定值，此外邊界粒子作用區(qū)域不完整，計算時存在由“表面效應”引發(fā)的誤差，同時計算時沒有考慮非局部作用程度與作用距離之間的關系導致定量計算時精度受到影響。

本文結合前人的基礎，將近場動力學中的微模量與彈性常數(shù)建立聯(lián)系，通過剛度等效的方式建立線性方程組，并且利用求解不定方程組最小二乘最小范數(shù)的方式來矯正微模量，在此基礎上引入可以反映非局部作用程度的核函數(shù)來提高計算精度，利用優(yōu)化后的方法對二維脆性材料在荷載作用下的變形及裂紋擴展過程進行了模擬。結果表明，本文提出的方法可以有效地降低邊界處的誤差并減少誤差產(chǎn)生的范圍，在模擬裂紋擴展時會有更好的收斂效果，較好的應用于固體材料的變形與破壞。

1 理論與假定

1.1 近場動力學理論

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

1.2 微彈性模型

1.3 小變形假定

圖2 小變形假設Fig. 2 Small deformation assumption

可發(fā)現(xiàn)模量與文獻[26]中相同，這可以證明本文小變形假定與鍵基近場動力學理論的一致性。c僅由體積模量K決定，導致該模型適用具有局限，另外在邊界區(qū)域由于積分區(qū)域不完整會產(chǎn)生較大的計算誤差。

2 非局部最小二乘優(yōu)化

2.1 非局部優(yōu)化

已證明模量c在傳統(tǒng)模型中為常量，通常非局部作用程度會隨距離而改變，使用常量c會導致定量計算時精度受到影響，Huang 等[28]針對該問題引入反映非局部作用程度的核函數(shù)從而提高定量計算的精度。本節(jié)在前人基礎上，在小變形假設推導的理論基礎上引入非局部核函數(shù)，以提高后續(xù)LSM 優(yōu)化計算精度，式(9)作用力中的可重新表示為：

2.2 最小二乘方法優(yōu)化

本節(jié)在小變形假定及非局部優(yōu)化的基礎上，將近場動力學模量與連續(xù)介質力學中彈性常數(shù)建立聯(lián)系，并利用最小二乘方法優(yōu)化鍵常數(shù)以此減少邊界誤差甚至突破適用局限性。式(12)可以表示為：

在經(jīng)典連續(xù)介質力學二維問題中，應力-應變及剛度張量的一般形式為：

利用式(42)可建立傳統(tǒng)連續(xù)介質力學剛度矩陣與近場動力學鍵模量之間的關系從而描述材料的變形特性，為后續(xù)裂紋擴展模擬奠定基礎，類似式(43)可以將所有模量之間的關系表達為：

最小范數(shù)最小二乘解可采用Moore-Penrose 偽逆矩陣求解：

上述問題可采用二次規(guī)劃(QR)處理，二次規(guī)劃是非線性規(guī)劃中的特殊規(guī)劃問題，在約束最小二乘問題求解方面效果較好。以拉格朗日方法為例求解等式約束二次規(guī)劃問題，考慮如下問題：

3 數(shù)值驗證

3.1 單軸荷載

為驗證本文方法的有效性，本節(jié)對二維平板在單軸荷載條件下的變形進行模擬，并將結果與理論及經(jīng)典矯正結果做出對比。矩形板的尺寸及所受荷載如圖3 所示。

圖3 單軸加載模型 /mm Fig. 3 Uniaxial loading model

試件尺寸1 m×0.5 m，采用平面應力方式計算，試件彈性模量取E=20 GPa ，泊松比取ν=1/3，密度 取 ρ=2800 kg/m3，拉 伸 荷 載 取p0=20 MPa，以中心為坐標原點，試件在荷載作用下x及y方向位移的理論解為：

將模型離散為 200×100節(jié)點，節(jié)點間距為Δx=0.005 m ， 半徑取 δ=3.015Δx，不采用任何矯正措施的結果與理論解的誤差百分比如圖4 所示。

圖4 無矯正計算結果相對誤差Fig. 4 Relative error of uncorrected calculation results

由圖4 可知不利用任何矯正措施的模型存在較大的誤差，特別是在角落及邊界區(qū)域，在加載方向及非加載方向最大誤差分別為21.33%與 32.05%，誤差主要集中在角落處及加載端處。

邊界產(chǎn)生大量誤差的主要原因是邊界附近的節(jié)點鄰域不完整。目前較為常用的方法是利用應變能密度等價引入修正系數(shù)，利用該方法計算結果與誤差百分比分別如圖5 與圖6 所示。

圖5 應變能密度計算結果Fig. 5 Calculation results by strain energy density

圖6 應變能密度相對誤差Fig. 6 Relative error by strain energy density

計算解與理論解的誤差百分比為了方便對比誤差百分比低于0.1%時不考慮。由圖6 可知，修正之后的計算結果較大改善，在角落與加載端部的誤差顯著降低，中部誤差基本消失，加載方向的誤差最大為1.72%，在非加載方向的誤差最大為12.73%。

而利用本文提出的方法計算結果與誤差百分比分別如圖7 與圖8 所示。不同方法計算結果的相對誤差對比如表1 所示。

表1 最大相對誤差對比表Table 1 Strain softening stress-strain curve of rock

圖7 LSM 單軸加載結果Fig. 7 Uniaxial calculation results by LSM

圖8 LSM 單軸相對誤差Fig. 8 Uniaxial relative error by LSM

本文方法計算的相對誤差與經(jīng)典方法相比，在加載方向上最大相對誤差沒有較大的改善，但是誤差產(chǎn)生的范圍所減小。此外在非加載方向相比應變能密度方法有較大的改善，從12.73%降至4.35%并且產(chǎn)生誤差的范圍也有所降低，說明本文提出方法在模擬結構變形方面的可行性。

3.2 雙軸荷載

為進一步驗證方法的有效性，雙軸荷載條件下的變形進行模擬，尺寸及所受荷載如圖9 所示。

圖9 雙軸加載模型 /mmFig. 9 Biaxial loading model

兩側拉伸荷載取px=20 MPa，上下拉伸荷載取py=15 MPa，位移理論解可用疊加法得：

模型離散及節(jié)點間距等參數(shù)與3.1 節(jié)一致，經(jīng)典應變能密度等價計算誤差如圖10 所示。

圖10 雙軸應變能密度相對誤差Fig. 10 Biaxial relative error by strain energy density

經(jīng)典方法在兩個加載方向的誤差最大分別為1.57%與11.71%，對比圖6 與圖10，雙軸加載條件下邊界的最大誤差相比單軸加載有所降低產(chǎn)生誤差的范圍有較小的增加。利用本文的方法模擬結果與誤差分布如圖11 與圖12 所示。

圖11 LSM 雙軸加載結果Fig. 11 Biaxial calculation results by LSM

本文方法在兩個加載方向的誤差最大分別為1.53%與3.29%。對比圖8 與圖12，雙軸最大相對誤差相比單軸有所降低，并且對比圖10 與圖12，最大誤差比經(jīng)典方法具有改善，尤其在最小主應力方向，且產(chǎn)生誤差的范圍也有所減少，進一步驗證了本文方法的有效性。

圖12 LSM 雙軸相對誤差Fig. 12 Biaxial relative error by LSM

3.3 裂紋擴展

為驗證本文方法研究材料斷裂的有效性，本節(jié)對含裂紋的脆性材料在荷載作用下的動態(tài)破裂過程進行模擬并將結果與文獻[29]對比，尺寸及所受荷載如圖13 所示。

圖13 裂紋擴展模型[29] /mmFig. 13 Crack growth model[29]

參數(shù)與文獻[29]相同取彈性模量E=72 GPa，泊松比ν=0.22 ，密度 ρ =2440 kg/m3，兩側拉伸荷載p0=14 MPa ，斷裂能G=135 J/m2。

為探究相同計算成本下的收斂結果，模型離散選取與一組參考文獻類似的節(jié)點數(shù) 200×80，間距為 Δx=0.5 mm ， 半徑取 δ=3.015Δx，不同時步nt下裂紋擴展情況以及裂紋擴展速度隨時間的變化對比分別如圖14 與圖15 所示。

圖14 不同時步損傷分布Fig. 14 Damage distribution at different time steps

圖15 裂紋擴展速度對比Fig. 15 Crack growth speed comparison

對比本文模擬結果與文獻[29]結果可知本文優(yōu)化的方法可以較好的反映脆性材料的裂紋擴展過程，捕捉到裂紋分叉的現(xiàn)象。此外由圖15 可知，在大致相等的計算成本下本文方法得到的裂紋擴展速度收斂效果相比文獻[29]有所提升，這其中表現(xiàn)在收斂速度與收斂結果上，利用本文方法在同等計算成本下可更接近實驗測得值。

4 結論

本文建立了基于非局部最小二乘優(yōu)化的近場動力學數(shù)值模型并對脆性材料的變形及破壞進行了模擬，取得以下結論：

(1) 基于最小二乘最小范數(shù)最優(yōu)解的近場動力學方法可以較好模擬固體材料的變形與破壞。與經(jīng)典的近場動力學方法相比，本文方法可以較好的降低邊界處的最大相對誤差及產(chǎn)生誤差的范圍，同時可以較為準確的反映脆性材料的裂紋擴展過程。

(2) 在荷載條件下，利用近場動力學方法所模擬的結果與理論結果所產(chǎn)生的誤差在兩個方向分別存在于角落與最大主應力邊界，并且雙軸加載條件下產(chǎn)生的最大誤差相比單軸條件有所降低。