基于非局部微分算子的近場(chǎng)動(dòng)力學(xué)及固體材料應(yīng)力數(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ái)，國(guó)內(nèi)外學(xué)者針對(duì)此問(wèn)題已取得了豐碩研究成果[1-4]，采用理論方法探究材料破壞機(jī)制時(shí)通常需要一定簡(jiǎn)化，無(wú)法很好地反映真實(shí)的情況，利用試驗(yàn)手段存在成本較高、時(shí)間周期長(zhǎng)等困難。隨著計(jì)算機(jī)技術(shù)發(fā)展，數(shù)值模擬手段越來(lái)越多的被用來(lái)研究材料的力學(xué)特性。

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

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

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

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

1 基本理論

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

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

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

1.2 微彈性模型

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

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

2 非局部微分算子

2.1 積分方程

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

2.2 應(yīng)力-應(yīng)變表達(dá)

由近場(chǎng)動(dòng)力學(xué)與經(jīng)典Navier 平衡方程對(duì)應(yīng)關(guān)系[28]，這里將近場(chǎng)動(dòng)力學(xué)中某點(diǎn)應(yīng)力分量表示為：

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

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

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

3.1 平板拉伸

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

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

3.2 含圓孔平板

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

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

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

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

3.3 含裂隙平板

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

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

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

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

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

3.4 裂紋擴(kuò)展過(guò)程

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

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

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

4.1 離散間距影響

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

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

4.2 泰勒項(xiàng)數(shù)影響

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

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

4.3 權(quán)函數(shù)影響

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

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

5 結(jié)論

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

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

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