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

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

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

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

1 黏彈性近場(chǎng)動(dòng)力學(xué)的理論框架

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

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

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

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

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

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

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

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

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

在黏彈性近場(chǎng)動(dòng)力學(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?條件下的動(dòng)態(tài)極限鍵拉伸。

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

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

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

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

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

任何變形模式都可以分解為2種基本的變形模式：體積變形和純剪切變形。對(duì)于體積變形模式，只有正應(yīng)變?chǔ)旁诎逯挟a(chǎn)生，可得到體積變形模式下微彈性參數(shù)κ與微黏性參數(shù)τ。

同理，對(duì)于純剪切變形模式，當(dāng)物體變形時(shí)，只有工程剪應(yīng)變?chǔ)迷诎逯挟a(chǎn)生，可以得到純剪切變形模式下的微彈性參數(shù)κ與微黏性參數(shù)τ。

由于2個(gè)材料點(diǎn)之間的彈簧力和阻尼力始終沿其軸向，是一個(gè)一維分量，因此，在黏彈性近場(chǎng)動(dòng)力學(xué)中應(yīng)該只有一個(gè)微彈性參數(shù)和一個(gè)微黏性參數(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 語(yǔ)言編寫(xiě)了黏彈性近場(chǎng)動(dòng)力學(xué)程序，給出2個(gè)算例：動(dòng)荷載作用下裂紋的分叉模擬和含圓孔各向同性平板的拉伸破壞。

3.1 拉伸動(dòng)荷載作用下裂紋的分叉模擬

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

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

3.2 含圓孔平板的動(dòng)荷載拉伸

數(shù)值試驗(yàn)研究了在板中心具有圓孔的二維混凝土板在動(dòng)態(tài)拉伸荷載下的模擬，圓孔的直徑是10mm，施加在板上的幾何參數(shù)和荷載條件如圖5所示。黏彈性近場(chǎng)動(dòng)力學(xué)材料參數(shù)與3.1節(jié)算例的材料參數(shù)保持一致，材料點(diǎn)的邊長(zhǎng)和域半徑也與3.1節(jié)的算例相同。荷載通過(guò)位移邊界條件施加在板的兩端邊界區(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

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

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

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

4 結(jié)論

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

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

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

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

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

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

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

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