基于非局部LSM 優(yōu)化的近場(chǎng)動(dòng)力學(xué)及脆性材料變形模擬

馬鵬飛，李樹(shù)忱，王修偉，王曼靈

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

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

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

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

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

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

1 理論與假定

1.1 近場(chǎng)動(dòng)力學(xué)理論

圖1 物質(zhì)點(diǎn)相互作用Fig. 1 Interaction of material points

1.2 微彈性模型

1.3 小變形假定

圖2 小變形假設(shè)Fig. 2 Small deformation assumption

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

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

2.1 非局部?jī)?yōu)化

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

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

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

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

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

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

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

3 數(shù)值驗(yàn)證

3.1 單軸荷載

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

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

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

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

圖4 無(wú)矯正計(jì)算結(jié)果相對(duì)誤差Fig. 4 Relative error of uncorrected calculation results

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

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

圖5 應(yīng)變能密度計(jì)算結(jié)果Fig. 5 Calculation results by strain energy density

圖6 應(yīng)變能密度相對(duì)誤差Fig. 6 Relative error by strain energy density

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

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

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

圖7 LSM 單軸加載結(jié)果Fig. 7 Uniaxial calculation results by LSM

圖8 LSM 單軸相對(duì)誤差Fig. 8 Uniaxial relative error by LSM

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

3.2 雙軸荷載

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

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

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

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

圖10 雙軸應(yīng)變能密度相對(duì)誤差Fig. 10 Biaxial relative error by strain energy density

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

圖11 LSM 雙軸加載結(jié)果Fig. 11 Biaxial calculation results by LSM

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

圖12 LSM 雙軸相對(duì)誤差Fig. 12 Biaxial relative error by LSM

3.3 裂紋擴(kuò)展

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

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

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

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

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

圖15 裂紋擴(kuò)展速度對(duì)比Fig. 15 Crack growth speed comparison

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

4 結(jié)論

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

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

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