劉一鳴+黃丹+秦洪遠

摘要： 構(gòu)建考慮混凝土拉壓異性和宏觀斷裂特征的混凝土類材料非局部態(tài)型近場動力學(xué)本構(gòu)模型，并通過引入動態(tài)松弛、系統(tǒng)失衡判斷和力邊界等效等算法，構(gòu)建適于分析混凝土類材料和結(jié)構(gòu)變形破壞過程的態(tài)型近場動力學(xué)數(shù)值模擬體系.通過分組模擬和定量計算，分析算法的收斂性、計算精度和效率等問題；在此基礎(chǔ)上開展含不同角度中心裂紋混凝土板的破壞模擬.

關(guān)鍵詞： 混凝土； 裂紋擴展； 態(tài)型近場動力學(xué)； 非局部模型

中圖分類號： O346.1 文獻標志碼： A

0 引 言

材料和結(jié)構(gòu)的破壞機制及其數(shù)值模擬是計算力學(xué)研究的經(jīng)典難題，也是諸多工程領(lǐng)域關(guān)注的重點.傳統(tǒng)的有限元等數(shù)值方法由于在連續(xù)介質(zhì)框架下求解偏微分方程，所以在分析破壞問題時必須預(yù)先知道裂紋的存在與否及其位置和尺寸，計算時又需要判斷裂紋是否擴展及擴展路徑，并重新剖分網(wǎng)格，具有一定的復(fù)雜性.邊界元、擴展有限元、非連續(xù)有限元等方法及虛擬裂紋閉合等措施可以較好地處理不連續(xù)問題，但由于依然在連續(xù)性假設(shè)框架下進行不連續(xù)區(qū)域的特殊處理，模擬復(fù)雜的諸如三維群裂紋、多尺度動力破壞等問題時，依然有待進一步研究.

近年來，基于非局部積分思想的近場動力學(xué)（PeriDynamics，PD）方法[1-2]憑借其不需要求解空間微分方程而在模擬大變形及裂紋擴展、爆炸和沖擊破壞等強不連續(xù)力學(xué)問題方面的突出優(yōu)勢[3-5]，成為計算力學(xué)和相關(guān)領(lǐng)域研究的熱點.針對建筑工程中最廣泛使用的混凝土材料和結(jié)構(gòu)，SILLING[6]，DEMMIE[7]和KILIC等[8-9]采用近場動力學(xué)方法模擬簡單混凝土梁和柱的沖擊破壞、失穩(wěn)等問題，GERSTLE等[10-12]在模擬混凝土梁的基礎(chǔ)上，還進一步分析近場動力學(xué)模型參數(shù)對計算精度的影響.本課題組也曾根據(jù)近場動力學(xué)思想建模分析常規(guī)混凝土構(gòu)件的拉、壓和沖擊破壞過程.[13-15]

然而，已有的相關(guān)工作主要基于常規(guī)“鍵型”近場動力學(xué)模型，不論是采用早期的單參數(shù)微彈脆性模型[6-9]，還是應(yīng)用改進后的微極鍵型近場動力學(xué)模型[11-12]，均無法真正滿足實際混凝土材料和結(jié)構(gòu)破壞模擬特別是定量計算分析的需要（例如，固有的泊松比限制問題[1，8]，三維模擬時單參數(shù)微觀脆性模型中泊松比限定于0.25，而微極模型[11-12]的泊松比則必須小于0.25），對材料的泊松比具有一定的限制.為修正包括泊松比限制在內(nèi)的鍵型近場動力學(xué)模型的缺陷，SILLING等[16-17]進一步發(fā)展態(tài)型近場動力學(xué)理論，近年來已成為關(guān)注的熱點[18-19]，但其相關(guān)工作主要圍繞各種傳統(tǒng)本構(gòu)模型和態(tài)型近場動力學(xué)模型之間的轉(zhuǎn)化來開展，對于混凝土材料和結(jié)構(gòu)破壞這一工程實際問題的態(tài)型近場動力學(xué)模擬尚未見諸報道.

本文基于態(tài)型近場動力學(xué)理論，構(gòu)建考慮混凝土拉壓異性和宏觀斷裂特征的彈脆性態(tài)型近場動力學(xué)本構(gòu)模型，并且通過引入動態(tài)松弛、系統(tǒng)失衡判斷、力邊界等效等算法，構(gòu)建完整的統(tǒng)一求解混凝土材料和結(jié)構(gòu)變形破壞問題的數(shù)值體系.通過對混凝土板單軸拉伸問題的定量計算和分組模擬，分析方法的收斂性、計算精度和效率.在此基礎(chǔ)上，應(yīng)用本文的模型和方法模擬含不同角度中心裂紋混凝土板的裂紋擴展與破壞過程.

1 態(tài)型近場動力學(xué)理論和模型

1.1 態(tài)型近場動力學(xué)理論

將空間物質(zhì)視為由帶質(zhì)量、有代表性體積的系列物質(zhì)點x組成，其

僅與其近場范圍H（|ξ|<δ）內(nèi)的其他物質(zhì)點相互作用，見圖1.圖1中：x和x′分別表示參考構(gòu)形中2個物質(zhì)點位置矢量；u和u′分別表示當(dāng)前構(gòu)形中2個物質(zhì)點位移矢量；y=x+u和y′=x′+u′分別表示當(dāng)前構(gòu)形中2個物質(zhì)點的坐標矢量；ξ=x′-x表示2個物質(zhì)點的相對位置矢量；η=u′-u表示2個物質(zhì)點的相對位移矢量.

1.2 二維常規(guī)態(tài)型近場動力學(xué)本構(gòu)模型

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

2.1 局部阻尼引入

為便于應(yīng)用近場動力學(xué)方法求解準靜力問題，根據(jù)準靜力結(jié)果驗證模型的可靠型，引入經(jīng)典力學(xué)中的動態(tài)松弛法，在運動方程中引入局部阻尼項[8，15]，將物質(zhì)點運動方程轉(zhuǎn)化為

2.2 數(shù)值離散

對平面問題采用四邊形均勻離散（單位厚度），離散間距為Δx，則物質(zhì)點的代表面積為VΔx=Δx2，對于物質(zhì)點xi考慮近場范圍內(nèi)（xj-xi≤δ）的相互作用

2.3 數(shù)值算法

近場動力學(xué)方法的顯式動力學(xué)實現(xiàn)涉及空間離散方式、時間差分格式及數(shù)值積分方法，數(shù)值積分方法與空間離散方式相關(guān).本文采用中心差分對時間序列進行離散，

2.4 系統(tǒng)失衡力準則

3 數(shù)值算例

為驗證本文模型和算法，首先通過對混凝土板的單向拉伸變形進行定量計算分析，進而分別對含不同角度中心裂紋的混凝土板受拉伸時的裂紋擴展和破壞過程進行分析.根據(jù)實際情況，將問題簡化為平面應(yīng)力問題分析.

對粒子系統(tǒng)的力邊界條件施加問題是粒子類方法中的難點.對于長為L，寬為W的矩形板（見圖2），本文在近場動力學(xué)模擬時在模型的加載邊界上增加寬度為δ的加載區(qū)域，將外力載荷等效為加載區(qū)物質(zhì)點體力施加于加載區(qū)物質(zhì)點上，使模型邊界上應(yīng)力等效邊界條件一致.這樣做還可以消除往常近場動力學(xué)計算中常見的邊界弱化問題.

3.1 板拉伸定量分析

取200 mm×200 mm混凝土方板，密度ρ=2 400 kg/m3，彈性模量E=30 GPa，泊松比v=0.2，體積模量和剪切模量分別為K=16.67 GPa，G=12.5 GPa.在豎直方向施加均布拉伸載荷σ=3 MPa.模型離散情況見圖2，離散物質(zhì)點間距Δx=1 mm，模型共包含40 600個物質(zhì)點，其中邊界物質(zhì)點600個，近場范圍取δ取3Δx.使用控制變量法，通過改變時間步長、人工阻尼因數(shù)C等計算參數(shù)，觀察迭代過程中板下邊緣中點y方向位移的收斂變化情況.為方便觀察位移隨時間步變化情況，將收斂容差-設(shè)為0，程序根據(jù)最大迭代步終止計算.

（1）時間步長取Δt=2×10-7s，改變阻尼系數(shù)C時的位移收斂情況見圖3.y方向的位移斂于1.011 8×10-5 mm，與彈性力學(xué)解1×10-5 mm相比，誤差為1.18%.

由此可見，當(dāng)阻尼系數(shù)相對較小時，位移值在理論解附近劇烈振蕩且無法收斂；隨著阻尼系數(shù)增大，位移振蕩幅值降低，且最終收斂于理論解.當(dāng)阻尼系數(shù)達到1×108 kg/（m3·s）時，位移值呈單調(diào)收斂，表現(xiàn)出類似超阻尼現(xiàn)象.阻尼系數(shù)繼續(xù)增大，位移收斂逐漸變緩.合適地選取阻尼系數(shù)可以大幅提高準靜態(tài)求解的計算效率.根據(jù)本文分組比較與分析，建議將局部阻尼項阻尼系數(shù)取為5×107～5×108 kg/（m3·s）.

（2）取阻尼系數(shù)C=2×108 kg/（m3·s）不變，改變時間步長Δt進行分析，見圖4.隨著時間步長的逐漸增大，位移的收斂也加快，其中時間步長為3×10-7s時收斂速度最快，但當(dāng)時間步長取Δt=5×10-7s時，差分迭代出現(xiàn)不收斂.關(guān)于態(tài)型近場動力學(xué)計算中動力學(xué)迭代時間步長的理論推導(dǎo)目前還需進一步深入研究和討論.

通過對人工阻尼系數(shù)和時間步長尺寸的敏感性分析可知，在基于態(tài)型近場動力學(xué)的準靜力問題分析中，阻尼系數(shù)和時間步長的選取對計算效率影響非常關(guān)鍵.若選取合適的阻尼系數(shù)和時間步長，則計算將快速收斂（如圖3和4）.此外，依靠最大迭代數(shù)終止計算具有很強的不穩(wěn)定性且無法反映計算結(jié)果的收斂情況，因此引入收斂準則顯得極為重要.

根據(jù)式（19）的失衡力準則，分別選取不同收斂容差，進行分組計算.分別取C=1×108 kg/（m3·s）和Δt=2×10-7 s，計算結(jié)果見表1：當(dāng)收斂容差-=1×10-4時，計算誤差接近1%；當(dāng)收斂容差繼續(xù)減小時，計算誤差僅有細微變化，而收斂迭代步數(shù)不斷增加.綜合考慮計算的精度和效率，收斂容差一般取1×10-4即可.

3.2 含中心裂紋板的裂紋擴展模擬

在圖2所示的混凝土板中心位置設(shè)置20 mm×1 mm的初始裂紋，并在豎直方向施加均布拉伸載荷，計算混凝土材料的抗壓強度和抗拉強度.對含不同角度中心裂紋的混凝土板進行分級準靜力加載，觀察裂紋開展情況，并確定其開裂載荷.時間步長取為Δt=2×10-7s，局部阻尼系數(shù)取C=2×108 kg/（m3·s），收斂容差取為-=1×10-4.

在上述平面模型中，近場范圍為3倍物質(zhì)點間距即δ=3 mm，裂紋寬度小于近場范圍，在形成物質(zhì)點對時斷開所有跨越裂紋線的物質(zhì)點對形成初始裂紋.

對于典型的平面I-II復(fù)合型裂紋擴展問題，根據(jù)脆性材料的最大拉應(yīng)力準則可確定起裂角應(yīng)滿足sin θ0+（3cos θ0-1）cot β=0，其中θ0為起裂角，β為裂紋與加載方向夾角.當(dāng)β分別為30°，45°，60°和90°時，理論計算所得起裂角θ0分別為60°，53.1°，43.2°和0°.

逐級加載過程如下：在初始加載時，以0.030 MPa為載荷增量進行施加，每次確定系統(tǒng)平衡后施加下一增量步；在系統(tǒng)出現(xiàn)損傷后，退回一個增量步將載荷增量降低為0.003 MPa進行加載，模擬裂紋的穩(wěn)定擴展，直至最后失穩(wěn)擴展，系統(tǒng)無法平衡，不滿足式（19）.

圖5a～5d分別顯示含45°初始裂紋混凝土板中從裂尖出現(xiàn)損傷、演化成裂紋起裂直至最后裂紋失穩(wěn)擴展直至貫穿的全過程.根據(jù)分級加載算法可以確定，含45°初始裂紋混凝土板的最終破壞載荷為0.441 MPa，即載荷增加到0.441 MPa時，系統(tǒng)無法達到平衡狀態(tài)，無法承受下一載荷增量，隨著動力學(xué)方程的繼續(xù)迭代，裂紋自然的快速失穩(wěn)擴展.

其他含不同初始角度裂紋混凝土板的破壞過程與含45°初始裂紋板類似，均由彈性變形、損傷累積、起裂和穩(wěn)定擴展、失穩(wěn)擴展等階段組成，主要區(qū)別在于最終的破壞載荷以及裂紋擴展路徑，詳見圖6.由圖6可知，盡管初始裂紋的傾角不同，裂紋最終基本會沿著垂直于加載方向快速擴展，與試驗結(jié)果一致.各種情況下裂紋的起裂角也與根據(jù)最大周向拉應(yīng)力計算的理論結(jié)果較為吻合.此外，由模擬結(jié)果可知，當(dāng)初始裂紋與加載方向夾角β增大時，I型裂紋所占的比例增加，板的承載能力也隨之下降，最終破壞載荷減小.這一結(jié)果也符合實際情況.

4 結(jié) 論

（1）考慮混凝土類材料的拉壓異性和宏觀斷裂

特征，建立混凝土新的態(tài)型近場動力學(xué)本構(gòu)模型，并通過引入動態(tài)松弛、系統(tǒng)失衡判斷、力邊界等效等系列算法，構(gòu)建完整的適合于求解混凝土材料和結(jié)構(gòu)變形破壞過程的態(tài)型近場動力學(xué)數(shù)值體系.

（2）通過對典型混凝土板的準靜態(tài)定量變形分析，驗證本文模型和算法的收斂性、計算精度和計算效率，并討論局部阻尼、時間步長和收斂容差等計算變量的選取對態(tài)型近場動力學(xué)計算精度和效率的影響，給出合理取值范圍建議.

（3）通過分級加載實現(xiàn)混凝土類材料和結(jié)構(gòu)從彈性變形、損傷累積、裂紋穩(wěn)定擴展以及失穩(wěn)破壞全過程的連續(xù)模擬與定量分析，并成功開展含不同角度中心裂紋混凝土板的破壞過程模擬，模擬所得起裂角和裂紋擴展路徑與試驗結(jié)果和理論值吻合，對結(jié)構(gòu)破壞載荷的定量計算結(jié)果與變化規(guī)律符合實際情況.

參考文獻：

[1] SILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces[J]. Journal of the Mechanics and Physics of Solids， 2000， 48（1）：175-209. DOI： 10.1016/S0022-5096（99）00029-0.

[2] 黃丹， 章青， 喬丕忠， 等. 近場動力學(xué)方法及其應(yīng)用[J]. 力學(xué)進展， 2010， 40（4）： 448-459.

HUANG D， ZHANG Q， QIAO P Z， et al. A review on peridynamics（PD） method and its application[J]. Advance in Mechanics， 2010， 40（4）： 448-459.

[3] LAI X， REN B， FAN H， et al. Peridynamics simulations of geomaterial fragmentation by impulse loads[J]. International Journal for Numerical and Analytical Methods in Geomechanics， 2015， 39（12）：1304-1330. DOI： 10.1002/nag.2356.

[4] FAN H， BERGEL G L， LI S. A hybrid peridynamics-SPH simulation of soil fragmentation by blast loads of buried explosive[J].International Journal of Impact Engineering， 2016（87）： 14-27. DOI： 10.1016/j.ijimpeng.2015.08.006.

[5] SILLING S A， ASKARI A. Peridynamic model for fatigue cracking， SAND2014-18590[R]. Sandia Report. New Mexico， 2014.

[6] SILLING S A， DEMMIE P， WARREN T L. Peridynamic simulation of high-rate material failure， SAND2007-3464C[R]. Proceedings of 2007 Applied Mechanics and Materials conference. Austin TX， 2007.

[7] DEMMIE P N， SILLING S A. An approach to modeling extreme loading of structure using peridynamics[J]. Journal of Mechanics of Materials and Structures， 2007， 2（10）： 1921-1945.

[8] KILIC B. Peridynamic theory for progressive failure prediction in homogeneous and heterogeneous materials[D].Tucson： The University of Arizona， 2008.

[9] KILIC B， MADENCI E. Structure stability and failure analysis using peridynamic theory[J]. International of Journal of Non-linear Mechanics， 2009， 44（8）： 845-854.

[10] GERSTLE W， SAU N， SILLING S A. Peridynamic modeling of plain and reinforced concrete structures[C]// Proceedings of 18th International Conference Structural Mechanics in Reactor Technology（SMiRT18）. Beijing， 2005.

[11] GERSTLE W， SAU N， AGUILERA E. Micropolar peridynamic constitutive model for concrete[C]// Proceedings of 19th International Conference Structural Mechanics in Reactor Technology （SMiRT19）， Toronto，2007： 12-17.

[12] GERSTLE W H， SAU N， SAKHAVAND N. On peridynamic computational simulation of concrete structures[J]. ACI Special Publication， 2009， 65： 245-264.

[13] 沈峰， 章青， 黃丹， 等. 基于近場動力學(xué)理論的混凝土軸拉破壞過程模擬[J]. 計算力學(xué)學(xué)報， 2013， 30（S1）： 79-83. DOI： 10.7511/jslx2013z017.

SHEN F， ZHANG Q， HUANG D， et al. Damage and failure process of concrete structure under uni-axialtension based on peridynamics modeling[J]. Chinese Journal of Computational Mechanics， 2013， 30（S1）： 79-83. DOI： 10.7511/jslx2013z017.

[14] HUANG D， ZHANG Q， QIAO P Z. Damage and progressive failure of concrete structures using non-local peridynamic modeling[J].Science China： Technological Science， 2011， 54（3）： 591-596. DOI： 10.1007/s11431-011-4306-3.

[15] HUANG D， LU G， QIAO P. An improved peridynamic approach for quasi-static elastic deformation and brittle fracture analysis[J]. International Journal of Mechanical Sciences， 2015（94/95）： 111-122. DOI： 10.1016/j.ijmecsci.2015.02.018.

[16] SILLING S A， EPTON M， WECKNER O， et al. Peridynamic states and constitutive modeling[J]. Journal of Elasticity， 2007， 88（2）： 151-184. DOI： 10.1007/s10659-007-9125-1.

[17] SILLING S A. Linearized theory of peridynamic states[J]. Journal of Elasticity， 2010， 99（1）： 85-111. DOI： 10.1007/s10659-009-9234-0.

[18] LE Q V， CHAN W K， SCHWARTZ J. A two-dimensional ordinary， state-based peridynamic model for linearly elastic solids[J]. International Journal for Numerical Methods in Engineering 2014， 98（8）： 547-561. DOI： 10.1002/nme.4642.

[19] WU C T. Kinematic constraints in the state-based peridynamics with mixed local/nonlocal gradient approximations[J]. Computational Mechanics， 2014， 54（5）： 1255-1267. DOI： 10.1007/s00466-014-1055-8.