混凝土板裂紋擴(kuò)展的態(tài)型近場(chǎng)動(dòng)力學(xué)模擬

劉一鳴+黃丹+秦洪遠(yuǎn)

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

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

中圖分類號(hào)： O346.1 文獻(xiàn)標(biāo)志碼： A

0 引 言

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

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

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

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

1 態(tài)型近場(chǎng)動(dòng)力學(xué)理論和模型

1.1 態(tài)型近場(chǎng)動(dòng)力學(xué)理論

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

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

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

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

2.1 局部阻尼引入

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

2.2 數(shù)值離散

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

2.3 數(shù)值算法

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

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

3 數(shù)值算例

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

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

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ì)點(diǎn)間距Δx=1 mm，模型共包含40 600個(gè)物質(zhì)點(diǎn)，其中邊界物質(zhì)點(diǎn)600個(gè)，近場(chǎng)范圍取δ取3Δx.使用控制變量法，通過改變時(shí)間步長(zhǎng)、人工阻尼因數(shù)C等計(jì)算參數(shù)，觀察迭代過程中板下邊緣中點(diǎn)y方向位移的收斂變化情況.為方便觀察位移隨時(shí)間步變化情況，將收斂容差-設(shè)為0，程序根據(jù)最大迭代步終止計(jì)算.

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

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

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

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

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

3.2 含中心裂紋板的裂紋擴(kuò)展模擬

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

參考文獻(xiàn)：

[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] 黃丹， 章青， 喬丕忠， 等. 近場(chǎng)動(dòng)力學(xué)方法及其應(yīng)用[J]. 力學(xué)進(jìn)展， 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] 沈峰， 章青， 黃丹， 等. 基于近場(chǎng)動(dòng)力學(xué)理論的混凝土軸拉破壞過程模擬[J]. 計(jì)算力學(xué)學(xué)報(bào)， 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.