路德春，宋志強(qiáng)，王國盛

（北京工業(yè)大學(xué)城市與工程安全減災(zāi)教育部重點(diǎn)實(shí)驗(yàn)室，北京 100124）

在實(shí)際工程中，混凝土結(jié)構(gòu)通常是帶裂縫工作的，并且破壞過程具有強(qiáng)非連續(xù)特性。傳統(tǒng)連續(xù)介質(zhì)力學(xué)方法為數(shù)值計(jì)算提供了豐富的理論基礎(chǔ)，基于連續(xù)介質(zhì)力學(xué)的部分?jǐn)?shù)值方法為混凝土非連續(xù)變形的模擬提供了方法和可能性，例如擴(kuò)展有限元方法［1］、邊界元法［2］、有限差分法［3］、無網(wǎng)格法［4-5］等，成功解決了混凝土結(jié)構(gòu)的一系列工程問題。然而，上述方法的控制方程仍然采用微分方程，在處理混凝土裂縫開展帶來的非連續(xù)問題時，仍然會導(dǎo)致對裂紋擴(kuò)展的不合理預(yù)測。Silling［6］提出的近場動力學(xué)方法將物體離散為具有體積和質(zhì)量的材料點(diǎn)。所有材料點(diǎn)自身的位移共同決定了物體的變形，材料點(diǎn)之間的相互作用是由材料點(diǎn)的相對位移產(chǎn)生，這種相互作用（也被稱為鍵）的極限狀態(tài)決定了材料的強(qiáng)度。描述材料點(diǎn)運(yùn)動的控制方程采用空間積分格式，避免了在裂紋位移不連續(xù)處計(jì)算空間導(dǎo)數(shù)的無效性，從而為模擬不連續(xù)變形提供了新的思路和手段。目前，近場動力學(xué)已經(jīng)成為研究材料裂紋形成、發(fā)展和分叉等非連續(xù)問題的重要方法［7-14］?；炷猎诩虞d過程中會出現(xiàn)開裂現(xiàn)象，其開裂模式、變形和強(qiáng)度特性在很大程度上取決于加載速率及其梯度［15-19］，這也表現(xiàn)為黏性效應(yīng)、熱活化效應(yīng)和慣性效應(yīng)［20］。近場動力學(xué)處理非連續(xù)變形問題因其優(yōu)勢，已經(jīng)被廣泛用于模擬混凝土的動態(tài)開裂過程［21-26］。然而，原始近場動力學(xué)［6］材料點(diǎn)的運(yùn)動方程僅包含與慣性效應(yīng)對應(yīng)的加速度項(xiàng)，缺少與黏性效應(yīng)對應(yīng)的速度項(xiàng)。從而帶來了對動態(tài)荷載作用下混凝土變形預(yù)測不合理的問題。Silling［27-28］發(fā)現(xiàn)，在原始近場動力學(xué)運(yùn)動方程中引入與阻尼相關(guān)的速度項(xiàng)可以解決波傳播的不穩(wěn)定性問題，相關(guān)的工作在最近的研究中取得了長足進(jìn)展［14，29-30］。由于阻尼力是一種非保守力，因此使用僅包含保守力所做功的虛功原理推導(dǎo)近場動力學(xué)運(yùn)動方程的方案不再適用。要合理模擬材料的變形，必須基于能量守恒原理引入非保守力做功，發(fā)展一種新的嚴(yán)格遵循熱力學(xué)原理的近場動力學(xué)方法。另一方面，原始近場動力學(xué)的斷裂準(zhǔn)則也無法合理反映混凝土強(qiáng)度的速率相關(guān)特性［22，31-32］。混凝土的強(qiáng)度本質(zhì)上是材料本構(gòu)層面的內(nèi)生特性，應(yīng)通過強(qiáng)度或斷裂準(zhǔn)則來表征［19，33-36］。因此，在鍵中引入率相關(guān)斷裂準(zhǔn)則也是十分必要的。

本文目的是發(fā)展包含非保守力做功的近場動力學(xué)方法?；诨炷恋淖冃螜C(jī)理，確定材料點(diǎn)的黏彈性相互作用模式，發(fā)展黏彈性近場動力學(xué)材料點(diǎn)的運(yùn)動方程和鍵的斷裂準(zhǔn)則。對于材料點(diǎn)的運(yùn)動，通過包含非保守力所做功的能量方程獲得控制方程。對于鍵的破壞，利用混凝土動態(tài)單軸S 強(qiáng)度準(zhǔn)則發(fā)展的率相關(guān)斷裂準(zhǔn)則來描述。并通過數(shù)值算例驗(yàn)證和應(yīng)用所提出的近場動力學(xué)方法。

1 黏彈性近場動力學(xué)的理論框架

近場動力學(xué)的基本思想是將連續(xù)材料離散為具有體積V和質(zhì)量m的有限個材料點(diǎn)，任意材料點(diǎn)i與一定范圍H（i）內(nèi)的其他材料點(diǎn)都具有相互作用。材料點(diǎn)i的域內(nèi)所有材料點(diǎn)相互作用的集合控制著點(diǎn)i的運(yùn)動，該運(yùn)動產(chǎn)生的慣性力與總的相互作用力保持平衡，其運(yùn)動方程為［6］

式中：k為微彈性模量；s為鍵拉伸，其等于鍵變形后長度與初始長度的差與初始長度的比值；μ為歷史函數(shù)。

現(xiàn)有研究表明［20］，混凝土的動態(tài)承載力包括慣性承載力和真實(shí)動態(tài)強(qiáng)度。真實(shí)動態(tài)強(qiáng)度受熱活化效應(yīng)和黏性效應(yīng)控制，這2 種效應(yīng)可以分別用彈簧和牛頓黏壺模型來表征。受混凝土動態(tài)物理機(jī)制的啟發(fā)，建立了鍵型近場動力學(xué)材料點(diǎn)黏彈性相互作用模型，如圖1 所示。模型中的阻尼力是一種非保守力，用于實(shí)現(xiàn)系統(tǒng)能量耗散的功能，所發(fā)展的模型反映了由非保守力所驅(qū)動的能量耗散。

圖1 近場動力學(xué)黏彈性相互作用模型Fig.1 Viscoelastic interaction model of peridynamics

1.1 黏彈性運(yùn)動方程

基于發(fā)展的黏彈性相互作用模型，通過包括非保守力做功的哈密爾頓原理［39-40］推導(dǎo)近場動力學(xué)黏彈性運(yùn)動方程，其表達(dá)式為

式中：δ為變分符號；t為時間；Wnc為非保守力所做的功；T為總動能；U為所有保守力產(chǎn)生的總勢能。T、U表達(dá)式如式（4）：

新建立的運(yùn)動方程既包括由彈簧力和體力組成的保守力項(xiàng)，也包括由阻尼力組成的非保守力項(xiàng)。當(dāng)非保守力項(xiàng)等于零時，黏彈性運(yùn)動方程式退化為原始近場動力學(xué)運(yùn)動方程，它是發(fā)展的運(yùn)動方程一個特例。

1.2 混凝土率相關(guān)斷裂準(zhǔn)則

在黏彈性近場動力學(xué)中，同樣采用應(yīng)變率相關(guān)的最大拉應(yīng)變準(zhǔn)則作為材料點(diǎn)的斷裂準(zhǔn)則，結(jié)合斷裂準(zhǔn)則的鍵力密度矢量可表示為

式中：FS為彈簧的彈力；FD為黏壺的阻尼力；μ被重新定義為

f(s0，s?)=s0d是斷裂準(zhǔn)則，s0d為鍵拉伸率為s?條件下的動態(tài)極限鍵拉伸。

對于黏彈性近場動力學(xué)方法的率相關(guān)斷裂準(zhǔn)則，通過定義動態(tài)極限鍵拉伸來描述極限鍵拉伸隨加載速率的增加，鍵拉伸的動態(tài)增長因子定義為：DIFt=s0d/s0。對于強(qiáng)度的DIFs，目前已經(jīng)提出了許多率相關(guān)公式［20，33，41-43］。在這些公式中，Lu等［20］提出的動態(tài)單軸S準(zhǔn)則能夠合理描述混凝土的真實(shí)動態(tài)強(qiáng)度。通過類比S 準(zhǔn)則給出了DIFt的率相關(guān)公式，如式（12）：

式中：?為對數(shù)鍵拉伸率，?=lgs?；λ和χmax是2 個參數(shù)，其確定方法與Lu等［20］提出的方法相同。

2 材料參數(shù)的標(biāo)定方法

發(fā)展的黏彈性近場動力學(xué)中包含2 類材料參數(shù)：微彈性參數(shù)和微黏性參數(shù)，通過將黏彈性近場動力學(xué)和連續(xù)介質(zhì)力學(xué)之間的能量密度進(jìn)行等效，即可獲得近場動力學(xué)的材料參數(shù)。

在二維條件下，連續(xù)介質(zhì)力學(xué)的黏彈性本構(gòu)關(guān)系被用來計(jì)算耗散能密度，其方程如式（13）：

任何變形模式都可以分解為2種基本的變形模式：體積變形和純剪切變形。對于體積變形模式，只有正應(yīng)變ε在板中產(chǎn)生，可得到體積變形模式下微彈性參數(shù)κ與微黏性參數(shù)τ。

同理，對于純剪切變形模式，當(dāng)物體變形時，只有工程剪應(yīng)變γ在板中產(chǎn)生，可以得到純剪切變形模式下的微彈性參數(shù)κ與微黏性參數(shù)τ。

由于2個材料點(diǎn)之間的彈簧力和阻尼力始終沿其軸向，是一個一維分量，因此，在黏彈性近場動力學(xué)中應(yīng)該只有一個微彈性參數(shù)和一個微黏性參數(shù)。聯(lián)立式（19）和式（20），得到K=2G、η1=2η2。根據(jù)K和G以及η1和η2之間的轉(zhuǎn)換關(guān)系，可以得出彈性與黏性泊松比在二維情況下始終都為常數(shù)1/3。因此，微彈性參數(shù)和微黏性參數(shù)可以統(tǒng)一表示為

3 數(shù)值算例

采用FORTRAN 90 語言編寫了黏彈性近場動力學(xué)程序，給出2個算例：動荷載作用下裂紋的分叉模擬和含圓孔各向同性平板的拉伸破壞。

3.1 拉伸動荷載作用下裂紋的分叉模擬

圖2為帶初始裂縫的混凝土平板在不同沖擊應(yīng)力條件下的開裂模式，對應(yīng)的沖擊應(yīng)力大小分別為σ1= 0.2 MPa，σ2= 0.5 MPa，σ3= 0.6 MPa，σ4=1.4 MPa［44］。與試驗(yàn)相同工況的數(shù)值模型如圖3 所示，離散材料點(diǎn)的直徑d=0.25 mm，每個材料點(diǎn)的域半徑選擇為Δ=3d=0.75mm。混凝土的材料參數(shù)為楊氏模量E=19.2GPa，黏性參數(shù)η=8.28 x 10-22N·s·m-2，ρ= 2 700 kg·m-3，S 準(zhǔn)則的參數(shù)λ=1.4，χmax=4。邊界的體積校正和表面效應(yīng)的處理與Silling等［37］和Madenci等［38］的處理方法相同。荷載通過應(yīng)力邊界條件施加在板的兩端邊界區(qū)域上，荷載條件與試驗(yàn)相同。

圖2 試驗(yàn)中不同應(yīng)力狀態(tài)下板的破壞模式Fig.2 Failure modes of plates in different stress states in experiment

圖3 模型幾何尺寸及荷載條件Fig.3 Geometric dimensions and load conditions of the model

在σ2=0.5MPa 的應(yīng)力狀態(tài)下，黏彈性近場動力學(xué)模擬的裂紋擴(kuò)展過程材料點(diǎn)的損傷如圖4b 所示，在加載的初期，裂紋從初始裂紋的尖端起裂并且新裂紋沿著初始裂紋方向擴(kuò)展。當(dāng)擴(kuò)展達(dá)到一定長度之后，裂紋出現(xiàn)了分叉現(xiàn)象，由1 條裂紋變成了2條裂紋，2 條分叉裂紋分別沿著與初始裂紋呈大約30°的方向向右邊界擴(kuò)展，最終達(dá)到右邊界，整個試件破壞。模擬的結(jié)果與其他學(xué)者模擬結(jié)果［45］吻合較好，進(jìn)而驗(yàn)證了發(fā)展的黏彈性近場動力學(xué)方法的有效性。在不同應(yīng)力條件下，裂紋的擴(kuò)展過程有很大區(qū)別。當(dāng)σ1=0.2MPa 時，混凝土板的新裂紋完全沿著初始裂紋的方向發(fā)展，近似一條直線并且沒有分叉，這種情況下的破壞模式與準(zhǔn)靜態(tài)破壞相同。當(dāng)應(yīng)力增加σ2=0.5MPa 時，出現(xiàn)了一次分叉，有2條主裂紋。當(dāng)應(yīng)力增加σ3= 0.6MPa 時，出現(xiàn)了2次分叉，有3條主裂紋。其中一次分叉形成的2條裂紋中有1 條產(chǎn)生二次分叉。當(dāng)應(yīng)力增加σ4=1.4MPa時，出現(xiàn)了3次分叉，產(chǎn)生4條主裂紋以及更多的小裂紋?；炷潦且环N率敏感的準(zhǔn)脆性材料，動態(tài)荷載會影響其裂紋的擴(kuò)展模式，并且隨著加載速率的增加，裂紋的分叉次數(shù)會逐漸增多，這是因?yàn)榧虞d速率越大，平衡輸入的能量所需要的裂紋區(qū)域更多。黏彈性近場動力學(xué)方法能很好地捕獲混凝土的率相關(guān)行為，并且對混凝土開裂等非連續(xù)問題能進(jìn)行合理描述。

圖4 不同應(yīng)力狀態(tài)下板的破壞模式Fig.4 Failure modes of plates in different stress states

3.2 含圓孔平板的動荷載拉伸

數(shù)值試驗(yàn)研究了在板中心具有圓孔的二維混凝土板在動態(tài)拉伸荷載下的模擬，圓孔的直徑是10mm，施加在板上的幾何參數(shù)和荷載條件如圖5所示。黏彈性近場動力學(xué)材料參數(shù)與3.1節(jié)算例的材料參數(shù)保持一致，材料點(diǎn)的邊長和域半徑也與3.1節(jié)的算例相同。荷載通過位移邊界條件施加在板的兩端邊界區(qū)域上，3 種加載速率分別為：v1= 5 m·s-1，v2=20 m·s-1，v3=50 m·s-1。在圖5中邊界的體積校正和表面效應(yīng)的處理與Silling 等［37］和Madenci等［38］的處理方法相同。

圖5 模型幾何尺寸及荷載條件Fig.5 Geometric dimensions and load conditions of the model

原始近場動力學(xué)和黏彈性近場動力學(xué)模擬的混凝土板在不同加載速率下的荷載-位移曲線與損傷破壞模式分別如圖6 和圖7 所示。對于黏彈性近場動力學(xué)，荷載-位移曲線的剛度受加載速率的影響十分明顯，隨著加載速率的增加而增加，并且極限承載能力也隨著加載速率而提高。不同加載速率下?lián)p傷范圍的區(qū)別也十分顯著，隨著加載速率增加，損傷破壞區(qū)域由一條非常細(xì)的直線逐漸變寬直至變成V形。這種現(xiàn)象的產(chǎn)生是因?yàn)榛炷林猩婕白冃蔚牟牧宵c(diǎn)的數(shù)量隨著加載速率的增加而增加，更多的材料點(diǎn)參與能量耗散過程。然而在不同加載速率下，原始近場動力學(xué)模擬得到的荷載-位移曲線的斜率基本沒有變化，極限承載力有部分提升，并不顯著，同時損傷破壞區(qū)域受加載速率影響的現(xiàn)象也十分微弱。通過對比得出，黏彈性近場動力學(xué)能夠合理捕獲混凝土剛度和強(qiáng)度都受加載速率影響的特性，并且能夠反應(yīng)混凝土非連續(xù)破壞的加載速率相關(guān)現(xiàn)象。

圖6 不同加載速度下原始近場動力學(xué)（PD）與黏彈性近場動力學(xué)（VPD）的荷載-位移曲線對比Fig.6 Comparison of load-displacement curves between peridynamics and viscoelastic peridynamics at different loading velocities

圖7 原始近場動力學(xué)與黏彈性近場動力學(xué)破壞模式對比Fig.7 Comparison of failure modes between original peridynamics and viscoelastic peridynamics

4 結(jié)論

在鍵型近場動力學(xué)方法基礎(chǔ)上，發(fā)展了黏彈性近場動力學(xué)材料點(diǎn)相互作用模型，根據(jù)材料點(diǎn)相互作用模型推導(dǎo)了黏彈性鍵型近場動力學(xué)運(yùn)動方程，通過與連續(xù)介質(zhì)力學(xué)能量密度等效獲得了材料參數(shù)的確定方法。利用動荷載下混凝土板的數(shù)值算例，驗(yàn)證了所提黏彈性近場動力學(xué)方法的有效性。主要結(jié)論如下：

（1）基于混凝土的變形機(jī)制，建立了材料點(diǎn)的黏彈性相互作用模型，基于哈密爾頓原理，推導(dǎo)了黏彈性運(yùn)動方程。該方程合理反映了混凝土在變形過程中由非彈性變形所產(chǎn)生的能量耗散以及混凝土材料的動態(tài)強(qiáng)度與變形行為。通過近場動力學(xué)和連續(xù)介質(zhì)力學(xué)之間的能量密度等效的方法，建立了近場動力學(xué)微彈性和微黏性參數(shù)的標(biāo)定方法。該方法所得到的微彈性參數(shù)與原始近場動力學(xué)微彈性參數(shù)一致。

（2）利用黏彈性近場動力學(xué)方法模擬了不同加載速率下的混凝土構(gòu)件的破壞模式，驗(yàn)證了該方法在模擬動荷載下混凝土強(qiáng)度、變形和裂紋擴(kuò)展模式方面的合理性。通過模擬混凝土帶有圓孔板在動態(tài)拉伸荷載下的力學(xué)性能，得到了受加載速率影響的破壞模式。

所提的黏彈性近場動力學(xué)方法能夠合理反映混凝土在非彈性變形下的能量耗散以及不同加載速率和應(yīng)力狀態(tài)下的裂紋擴(kuò)展和分叉，為混凝土動態(tài)力學(xué)性能的描述提供了一種有效的數(shù)值方法和工具。對于混凝土在更加真實(shí)的三維情況下的力學(xué)性能和破壞狀態(tài)的模擬還有待進(jìn)一步研究，并且對于復(fù)雜應(yīng)力狀態(tài)下的力學(xué)性能研究也亟待發(fā)展。

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

路德春：項(xiàng)目負(fù)責(zé)人，提供研究平臺、構(gòu)思論文、指導(dǎo)模型構(gòu)建及分析數(shù)據(jù)、修改論文。

宋志強(qiáng)：構(gòu)建模型、分析數(shù)據(jù)、實(shí)現(xiàn)程序、呈現(xiàn)結(jié)果及撰寫論文。

王國盛：項(xiàng)目負(fù)責(zé)人，提供理論及創(chuàng)新思路、提供編程幫助、構(gòu)思與修改論文。