陳希卓，禹海濤，朱建波，劉建鋒

（1.同濟(jì)大學(xué)土木工程學(xué)院，上海 200092；2.四川大學(xué)深地科學(xué)與工程教育部重點(diǎn)實(shí)驗(yàn)室，四川成都 610065；3.同濟(jì)大學(xué)土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092；4.深圳大學(xué)土木與交通工程學(xué)院，廣東深圳 518061）

材料的裂紋擴(kuò)展和破壞問(wèn)題一直都是固體力學(xué)研究和工程領(lǐng)域面臨的關(guān)鍵難題[1］。目前學(xué)者們已經(jīng)提出了眾多的力學(xué)模型和數(shù)值方法來(lái)模擬固體材料斷裂過(guò)程，如擴(kuò)展有限元[2-3］、內(nèi)聚力模型[4-5］和邊界元方法[6］等。然而，這些方法都需要引入特定的附加函數(shù)和斷裂準(zhǔn)則來(lái)描述不連續(xù)行為，無(wú)法有效模擬裂紋的自發(fā)萌生和擴(kuò)展過(guò)程[7-8］。此外，由于傳統(tǒng)的數(shù)值方法大都是基于連續(xù)性假設(shè)，因此在處理斷裂等非連續(xù)行為時(shí)必須面對(duì)導(dǎo)數(shù)奇異性問(wèn)題，從本質(zhì)上難以準(zhǔn)確模擬破壞問(wèn)題。

為了解決這一難題，Silling[7］提出了近場(chǎng)動(dòng)力學(xué)（PD）理論，采用考慮非局部作用的積分模型代替?zhèn)鹘y(tǒng)理論的微分模型，有效避免了傳統(tǒng)連續(xù)介質(zhì)力學(xué)方法在處理斷裂行為時(shí)的導(dǎo)數(shù)奇異性問(wèn)題。該方法實(shí)現(xiàn)了由連續(xù)到不連續(xù)、微觀到宏觀力學(xué)作用的統(tǒng)一描述，尤其在分析裂紋擴(kuò)展等非連續(xù)力學(xué)問(wèn)題時(shí)具有顯著優(yōu)勢(shì)。近些年來(lái)，PD方法已逐漸成為計(jì)算力學(xué)及相關(guān)領(lǐng)域的關(guān)注熱點(diǎn)，其主要可以分為鍵基和態(tài)基2種模型[8］。對(duì)比態(tài)基近場(chǎng)動(dòng)力學(xué)模型，鍵基模型的計(jì)算量更小，同時(shí)不存在零能模態(tài)的問(wèn)題，因此廣泛應(yīng)用于巖石和混凝土等準(zhǔn)脆性材料的斷裂行為模擬。朱其志等[9］運(yùn)用近場(chǎng)動(dòng)力學(xué)方法對(duì)含有預(yù)制裂隙的巖石類材料試件的單軸壓縮試驗(yàn)進(jìn)行數(shù)值模擬，分析了不同裂隙傾角對(duì)裂紋擴(kuò)展模式的影響。Rabczuk和Ren[10］基于帶對(duì)偶-近場(chǎng)作用的近場(chǎng)動(dòng)力學(xué)（DHPD）模型實(shí)現(xiàn)了巖石內(nèi)部裂紋擴(kuò)展和準(zhǔn)脆性斷裂過(guò)程模擬。Wang等[11］發(fā)展了一種新型共軛鍵基模型，模擬了單軸壓縮條件下預(yù)先含有裂隙的巖石試樣中裂隙的分叉和合并行為，揭示了單軸壓縮條件下巖石的裂紋擴(kuò)展機(jī)制。Gerstle等[12］提出了微極近場(chǎng)動(dòng)力學(xué)（MPPD）模型，采用Euler-Bernoulli梁模型描述物質(zhì)點(diǎn)之間的相互作用，有效模擬了混凝土構(gòu)件失穩(wěn)和破壞過(guò)程。Diana等[13］運(yùn)用微極近場(chǎng)動(dòng)力學(xué)分析了砂巖在三點(diǎn)彎曲試驗(yàn)中的力學(xué)特性，有效捕捉到了試件的I型和混合型斷裂過(guò)程。

上述近場(chǎng)動(dòng)力學(xué)模型在準(zhǔn)脆性材料斷裂過(guò)程的模擬中均取得了不錯(cuò)的效果，但其主要還是集中于二維斷裂問(wèn)題。實(shí)際上，準(zhǔn)脆性材料在三維條件下會(huì)出現(xiàn)更復(fù)雜的體破壞現(xiàn)象，裂紋的空間擴(kuò)展模式也會(huì)呈現(xiàn)較大差別[14］。而傳統(tǒng)的微極近場(chǎng)動(dòng)力學(xué)模型對(duì)復(fù)雜條件下材料三維斷裂行為的模擬精度不高，主要原因在于其無(wú)法保證不均勻應(yīng)變場(chǎng)下近場(chǎng)動(dòng)力學(xué)和傳統(tǒng)連續(xù)介質(zhì)力學(xué)的能量一致性以及缺乏相應(yīng)的三維斷裂準(zhǔn)則。

本文提出一種廣義微極近場(chǎng)動(dòng)力學(xué)（GMPD）模型，在傳統(tǒng)微極近場(chǎng)動(dòng)力學(xué)（MPPD）模型基礎(chǔ)上進(jìn)一步考慮三維狀態(tài)下鍵的軸向變形、切向變形、相對(duì)轉(zhuǎn)角之間的耦聯(lián)作用，可有效提高三維斷裂問(wèn)題的模擬精度。通過(guò)引入對(duì)應(yīng)于鍵拉、剪、彎力學(xué)行為的3種近場(chǎng)動(dòng)力學(xué)參數(shù)，以保證復(fù)雜荷載條件下近場(chǎng)動(dòng)力學(xué)與傳統(tǒng)連續(xù)介質(zhì)力學(xué)應(yīng)變能的一致性。另外建立基于能量的新型斷裂準(zhǔn)則，旨在實(shí)現(xiàn)準(zhǔn)脆性材料的三維斷裂過(guò)程模擬與準(zhǔn)確預(yù)測(cè)。

1 廣義微極近場(chǎng)動(dòng)力學(xué)模型

基于非局部作用的思想，近場(chǎng)動(dòng)力學(xué)理論假定物體內(nèi)的任一物質(zhì)點(diǎn)位置矢量xi與其周圍一定區(qū)域Hx內(nèi)的其他任意物質(zhì)點(diǎn)位置矢量xj之間存在相互作用，這種作用可以理解為鍵，而力通過(guò)鍵在物質(zhì)點(diǎn)間進(jìn)行傳遞。為了消除泊松比的限制，微極近場(chǎng)動(dòng)力學(xué)模型引入了Euler-Bernoulli梁模型來(lái)描述物質(zhì)點(diǎn)之間的相互作用，但是該模型忽略了鍵的剪切變形和轉(zhuǎn)動(dòng)慣量，無(wú)法對(duì)復(fù)雜加載條件下的力學(xué)行為尤其是三維斷裂問(wèn)題進(jìn)行準(zhǔn)確模擬。鑒于此，提出了一種廣義微極近場(chǎng)動(dòng)力學(xué)模型，引入了Timoshenko梁來(lái)模擬物質(zhì)點(diǎn)間的相互作用，進(jìn)一步提高模型的模擬精度與適用性。

1.1 控制方程

基于Timoshenko梁理論，xi在t時(shí)刻的運(yùn)動(dòng)控制方程可表示為

式中：ρ為物質(zhì)密度；u（xi，t）為xi的位移向量；Hx為xi的作用域；f為xi和xj間的相互作用力；θi和θj分別為xi和xj的相對(duì)轉(zhuǎn)角；b（xi，t）為作用于xi上的外部體力；I為轉(zhuǎn)動(dòng)慣量；A為鍵的截面面積；m為xi和xj間的彎矩相互作用；n（xi，t）為作用于xi上的外部彎矩；η和ξ分別為物質(zhì)點(diǎn)間的相對(duì)位置向量和相對(duì)位移向量，可表示為

基于局部坐標(biāo)系，鍵的對(duì)偶力函數(shù)f與對(duì)偶力矩m可分解為

式中：e1、e2和e3分別為局部坐標(biāo)系下沿x、y和z軸的單位向量。

引入分別對(duì)應(yīng)拉伸、剪切和彎曲剛度的鍵參數(shù)CN、Cθ和CM，并結(jié)合Timoshenko梁?jiǎn)卧芰Ψ匠?，可得到單根鍵上的力和彎矩在局部三維坐標(biāo)系下的表達(dá)為

1.2 應(yīng)變能與鍵參數(shù)

通過(guò)對(duì)鍵的微勢(shì)能w進(jìn)行域內(nèi)的積分，可以得到該模型的應(yīng)變能密度WPD。

式中：d為局部坐標(biāo)系下鍵的變形分量的統(tǒng)一表達(dá)；K為鍵的剛度矩陣。

建立如圖1所示的球面坐標(biāo)系OX′Y′Z′，三維條件下GMPD模型的應(yīng)變能密度可以進(jìn)一步表示為

圖1 球坐標(biāo)系下鍵的變形構(gòu)型Fig.1 A bond with the deformed configuration in global spherical coordinate system

式中：δ為作用域的半徑大??；α為鍵與球坐標(biāo)系中z′軸的夾角；β為鍵在Ox′y′平面的投影與x′軸的夾角。

基于近場(chǎng)動(dòng)力學(xué)應(yīng)變能WPD和傳統(tǒng)連續(xù)介質(zhì)力學(xué)應(yīng)變能WCM的等價(jià)關(guān)系可以得到三維條件下鍵參數(shù)的表達(dá)為

將式（5）和式（6）代入式（8），并通過(guò)局部到整體的坐標(biāo)變換可以得到GMPD模型的應(yīng)變能密度表達(dá)，如式（9）：

而傳統(tǒng)連續(xù)介質(zhì)力學(xué)的應(yīng)變能密度WCM為

1.3 斷裂準(zhǔn)則

為了有效捕捉準(zhǔn)脆性材料的破壞過(guò)程，提出了一種基于能量的新型鍵基斷裂準(zhǔn)則。在該模型中，假定當(dāng)穿越某一平面的鍵AB全部斷開(kāi)時(shí)，物質(zhì)點(diǎn)A的域內(nèi)形成相應(yīng)的宏觀裂紋，如圖2所示，因此材料的平面斷裂能應(yīng)等于所有斷開(kāi)鍵的應(yīng)變能總和，如式（13）：

圖2 斷裂演化過(guò)程示意Fig.2 Evaluation of fracture energy

展開(kāi)求解式（13），可以得到鍵的軸向變形、剪切變形以及相對(duì)轉(zhuǎn)角的限值s0、γ0和θ0，分別為

當(dāng)鍵的變形達(dá)到相應(yīng)的峰值后會(huì)自動(dòng)斷開(kāi)，即鍵兩端的物質(zhì)點(diǎn)不再產(chǎn)生相互作用，從而在宏觀上形成裂紋。這里引入標(biāo)量函數(shù)μ（xi，xj，t）描述物質(zhì)點(diǎn)間鍵是否發(fā)生破壞，如式（15）：

基于上述標(biāo)量函數(shù)，xi的損傷度d（xi，t）可以表示為

損傷度為0，代表物質(zhì)點(diǎn)完全無(wú)損傷；損傷度為1，表明物質(zhì)點(diǎn)周圍的鍵完全斷裂。

2 數(shù)值實(shí)現(xiàn)

為了求解GMPD模型中的控制方程，首先在空間上對(duì)求解域進(jìn)行離散化處理，在這里將求解域離散成物質(zhì)點(diǎn)，每個(gè)物質(zhì)點(diǎn)都占據(jù)一定空間體積且擁有一定物理性質(zhì)，離散后的運(yùn)動(dòng)控制方程可以表示為

式中：ΔVj為空間離散化處理后xi域內(nèi)的xj所占據(jù)的空間體積。

采用顯示中心差分方法求解上述控制方程，用有限差分近似代替位移對(duì)時(shí)間的導(dǎo)數(shù)，可以得到物質(zhì)點(diǎn)的加速度和速度的表達(dá)為

將式（19）代入式（17），可以得到位移的遞推求解格式為

將式（20）代入式（18），可以得到轉(zhuǎn)角的遞推求解格式為

該模型的邊界條件主要包括位移邊界和力邊界。其中外力P（x）是通過(guò)轉(zhuǎn)化為體力b（x）作用于最外層的物質(zhì)點(diǎn)Lr上，如圖3所示，其中Δ為物質(zhì)點(diǎn)間距。具體的轉(zhuǎn)換關(guān)系為

位移邊界主要是作用于物質(zhì)點(diǎn)外的虛擬邊界上，通過(guò)虛擬邊界層帶動(dòng)內(nèi)部物質(zhì)點(diǎn)一起運(yùn)動(dòng)。如圖3所示，虛擬邊界的尺寸和作用域半徑δ一致。

圖3 力和位移邊界條件Fig.3 Force and displacement boundary conditions

GMPD模型的完整數(shù)值求解流程如圖4所示。

圖4 數(shù)值求解流程Fig.4 Computational flowchart of proposed model

3 三維斷裂模擬

為了說(shuō)明GMPD模型對(duì)于三維斷裂問(wèn)題的適用性，基于該模型模擬單軸壓縮條件下準(zhǔn)脆性材料試件的準(zhǔn)靜態(tài)破壞過(guò)程，并通過(guò)與傳統(tǒng)微極模型以及文獻(xiàn)中已有試驗(yàn)結(jié)果進(jìn)行對(duì)比來(lái)驗(yàn)證模型的有效性。

選取含單條初始裂隙的立方體試件為對(duì)象，如圖5所示?；谖墨I(xiàn)[14］的試驗(yàn)參數(shù)，試件的彈性模量取E=3.81×109Pa，密度ρ=1.21g·cm-3，泊松比v=0.4，表面斷裂能G0=452N·m-1，初始裂紋半徑a=5mm，裂紋傾角α=30°，在兩側(cè)施加速率為2mm·min-1的軸向壓縮荷載。在近場(chǎng)動(dòng)力學(xué)模擬中，假定試件為均勻和各向同性的，采用間距Δ=0.001m的物質(zhì)點(diǎn)對(duì)試件進(jìn)行離散，作用域的大小取δ=3Δ。

圖5 壓縮荷載下含初始裂隙試件[14］（單位：mm）Fig.5 Single flawed samples under uniaxial compression[14](unit:mm)

圖6為加載過(guò)程中本模型預(yù)測(cè)的試件裂紋擴(kuò)展?fàn)顟B(tài)。由圖6a可知，預(yù)制裂紋邊緣最先開(kāi)始發(fā)生破壞并逐漸形成翼型裂紋。在翼型裂紋不斷延伸的同時(shí)，反翼型裂紋（方向和翼型裂紋相反的拉伸裂紋）也開(kāi)始在預(yù)制裂紋尖端萌生并擴(kuò)展（見(jiàn)圖6b、圖6c）。隨著荷載的逐漸增加，翼型裂紋邊緣產(chǎn)生了次生裂紋，裂紋沿水平方向延伸至試件邊緣造成材料破壞（見(jiàn)圖6d）。本文模擬所得裂紋類型及其擴(kuò)展路徑與試驗(yàn)結(jié)果[14］基本一致（如圖7），有效捕捉到了三維空 間內(nèi)產(chǎn)生的翼型裂紋、反翼型裂紋以及次生裂紋。

圖6 加載過(guò)程中GMPD模型預(yù)測(cè)的裂紋擴(kuò)展過(guò)程Fig.6 Predicted crack propagation in specimen under axial compression

圖7 壓縮荷載下試驗(yàn)觀測(cè)的試件裂紋擴(kuò)展過(guò)程[14］Fig.7 Experimental observation of crack development of specimen under axial compression[14]

為了進(jìn)一步說(shuō)明二維斷裂與三維斷裂模擬的差異性，圖8給出了相同荷載作用下傳統(tǒng)MPPD模型預(yù)測(cè)的二維裂紋的擴(kuò)展模式。

由圖8可見(jiàn)，二維裂紋發(fā)展主要以翼型裂紋為主，且擴(kuò)展方向與荷載施加方向基本一致。進(jìn)而，通過(guò)與圖6對(duì)比可以發(fā)現(xiàn)：三維條件下裂紋的擴(kuò)展過(guò)程更為復(fù)雜，會(huì)演化生成更多不同類型的裂紋，從而進(jìn)一步證明本文所建立的GMPD模型在模擬準(zhǔn)脆性材料三維裂紋擴(kuò)展問(wèn)題方面的準(zhǔn)確性和優(yōu)越性。

圖8 傳統(tǒng)MPPD模型預(yù)測(cè)的最終裂紋擴(kuò)展模式Fig.8 Final growth path obtained from the original MPPD model

4 結(jié)語(yǔ)

提出了一種廣義微極近場(chǎng)動(dòng)力學(xué)（GMPD）模型，可以有效模擬準(zhǔn)脆性材料的三維斷裂行為。模型采用Timoshenko梁來(lái)模擬物質(zhì)點(diǎn)間的相互作用，從而充分考慮鍵在三維受力條件下的軸向變形、切向變形、相對(duì)轉(zhuǎn)角以及三者之間的耦合作用。引入了可分別表征鍵的拉伸、剪切以及彎曲剛度的3個(gè)鍵參數(shù)，實(shí)現(xiàn)了任意變形場(chǎng)下近場(chǎng)動(dòng)力學(xué)和傳統(tǒng)連續(xù)介質(zhì)力學(xué)的能量一致。同時(shí)基于能量提出了一種新型斷裂準(zhǔn)則，給出了鍵軸向變形、切向變形以及相對(duì)轉(zhuǎn)角的臨界值，實(shí)現(xiàn)準(zhǔn)脆性材料的三維破壞模擬。

基于所提出的GMPD模型模擬了單軸壓縮荷載作用下試件的三維破壞過(guò)程，結(jié)果表明該模型可以有效捕捉復(fù)雜荷載條件下不同類型裂紋的萌生和擴(kuò)展，包括翼型裂紋、反翼型裂紋以及次生裂紋的發(fā)展過(guò)程。通過(guò)模型預(yù)測(cè)結(jié)果與試驗(yàn)結(jié)果對(duì)比，進(jìn)一步驗(yàn)證了模型的準(zhǔn)確性，表明本模型可以有效應(yīng)用于準(zhǔn)脆性材料的三維斷裂行為模擬，可為工程破壞問(wèn)題提供分析依據(jù)。

作者貢獻(xiàn)聲明：

陳希卓：模型構(gòu)建、數(shù)據(jù)分析呈現(xiàn)及論文撰寫(xiě)。

禹海濤：項(xiàng)目負(fù)責(zé)人、論文構(gòu)思、指導(dǎo)模型構(gòu)建及數(shù)據(jù)分析，論文修改。

朱建波：論文修改，提供試驗(yàn)數(shù)據(jù)。

劉建鋒：論文修改。