， ，志云(. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 604；. 大連海洋大學(xué) 海洋與土木工程學(xué)院，遼寧 大連 603)

混凝土材料通常指由水泥為膠凝材料，砂、石作骨料，與水(可含外加劑和摻合料)按一定比例配合，經(jīng)攪拌而得的水泥混凝土。在土木工程結(jié)構(gòu)中，混凝土材料得到了大量的應(yīng)用，其抗拉能力比抗壓能力弱很多，然而在一些極端情況下，混凝土材料可能會(huì)承受較大的拉力從而引起結(jié)構(gòu)破壞。巴西劈裂試驗(yàn)是測(cè)定巖石類材料抗拉強(qiáng)度的方法之一，也是國(guó)際巖石力學(xué)學(xué)會(huì)(international society for rock mechanics, ISRM)[1]和美國(guó)材料試驗(yàn)學(xué)會(huì)(american society for testing and materials, ASTM)[2]推薦的測(cè)量巖石材料抗壓強(qiáng)度的試驗(yàn)方法。由于混凝土與巖石力學(xué)性能相近，有許多學(xué)者利用巴西劈裂試驗(yàn)測(cè)量混凝土的抗拉強(qiáng)度[3-6]。Cundall和Strack于20世紀(jì)70年代基于牛頓第二定律提出了離散單元法(discrete element method, DEM)的概念。相比基于連續(xù)介質(zhì)力學(xué)的有限單元法，離散單元法對(duì)于模擬混凝土的開裂、裂紋擴(kuò)展等大變形有一定的優(yōu)勢(shì)。Qin等[7]基于質(zhì)點(diǎn)單元法，建立了不同應(yīng)變率下混凝土骨料、砂漿和界面三相細(xì)觀動(dòng)力學(xué)模型，研究表明應(yīng)變率越高，網(wǎng)狀裂紋產(chǎn)生的越多，試樣的力鏈產(chǎn)生更多分支，意味著在更高的應(yīng)變率下需要更多的能量達(dá)到失效。Ghazvinian等[8]用PFC2D模擬了含非貫通節(jié)理的巖石材料的剪切特性，離散元模擬表明宏觀剪切帶是張拉破壞產(chǎn)生的微裂紋擴(kuò)展引起的。Mechtcherine等[9]介紹了用離散單元法模擬混凝土流動(dòng)領(lǐng)域的發(fā)展，給出了參數(shù)估計(jì)的各種例子并與試驗(yàn)結(jié)果進(jìn)行了比較，顯示出離散單元法的潛力。Pan等[10]用PFC2D建立了混凝土三相細(xì)觀模型，進(jìn)行了單軸壓縮試驗(yàn)來(lái)研究堿集料反應(yīng)對(duì)混凝土力學(xué)性能的影響，結(jié)果表明這一模型可以有效預(yù)測(cè)堿集料反應(yīng)引起的混凝土內(nèi)部微裂紋的擴(kuò)展。Kazerani等[11]提出了一個(gè)檢驗(yàn)巖石強(qiáng)度和失效的離散元模型，兩個(gè)相鄰顆粒的邊界被視為柔性接觸，其本構(gòu)關(guān)系控制材料的斷裂性能，基于單軸/三軸壓縮和巴西劈裂試驗(yàn)提出一套細(xì)觀參數(shù)并對(duì)比試驗(yàn)結(jié)果驗(yàn)證了方法的有效性。

離散元單元法與有限單元法相比，最大的區(qū)別是不可能通過(guò)直接賦值的形式實(shí)現(xiàn)對(duì)材料的宏觀基本物理力學(xué)參數(shù)的確定，僅取決于顆粒的細(xì)觀本構(gòu)模型參數(shù)。離散元只可以直接賦值顆粒的幾何特征(顆粒大小和分布)和顆粒間接觸的細(xì)觀力學(xué)參數(shù)(細(xì)觀本構(gòu)模型和剛度、摩擦系數(shù)等接觸力學(xué)特征)，這些細(xì)觀參數(shù)的改變意味著材料的宏觀力學(xué)特性的顯著改變。而如何定量確定混凝土材料離散元模型的細(xì)觀參數(shù)是制約該方法進(jìn)一步深入研究和工程實(shí)際應(yīng)用的瓶頸問題之一?；炷凉橇狭浇?jīng)常是毫米量級(jí)的，在該尺度上直接測(cè)量其細(xì)觀參數(shù)是極其困難的，通用的方法是試錯(cuò)法(hit-and-miss)，即不斷調(diào)整細(xì)觀參數(shù)反復(fù)計(jì)算直到獲得一組與宏觀響應(yīng)吻合較好的細(xì)觀參數(shù)。然而這樣的方法往往存在著參數(shù)調(diào)整缺乏方向性和目的性等問題，造成效率低下。本研究提出一種基于巴西劈裂試驗(yàn)的混凝土細(xì)觀本構(gòu)模型參數(shù)估計(jì)方法，即先建立試驗(yàn)件離散元模型，給出一組細(xì)觀參數(shù)估計(jì)值，接下來(lái)分別改變每個(gè)細(xì)觀參數(shù)的取值，利用線性回歸確定每個(gè)細(xì)觀參數(shù)對(duì)宏觀變量的影響系數(shù)，進(jìn)而確定試樣的最大承載力及其對(duì)應(yīng)的最大位移與細(xì)觀參數(shù)間的關(guān)系，求解一個(gè)約束優(yōu)化問題，即可求得一組計(jì)算結(jié)果與實(shí)際試驗(yàn)曲線較為接近的細(xì)觀參數(shù)。該方法的主要優(yōu)勢(shì)是只使用一組試驗(yàn)數(shù)據(jù)并進(jìn)行少量PFC2D計(jì)算，即可得到一組滿足工程需求的細(xì)觀參數(shù)。

1 PFC2D中的細(xì)觀本構(gòu)模型

在某一時(shí)刻，PFC2D程序根據(jù)牛頓第二定律確定顆粒的加速度，并通過(guò)差分法得到顆粒在這一時(shí)刻的位移。而將力和位移聯(lián)系起來(lái)使計(jì)算能夠循環(huán)進(jìn)行下去的關(guān)鍵條件，就是細(xì)觀本構(gòu)模型。PFC2D中常用的細(xì)觀本構(gòu)模型主要有三種：一是接觸剛度模型，主要描述顆粒相互接觸運(yùn)動(dòng)時(shí)法向和切向的擠壓力；二是滑移模型，主要描述顆粒間存在相互滑動(dòng)趨勢(shì)時(shí)的力學(xué)行為；三是粘結(jié)模型，主要描述黏性顆粒的抗拉和抗剪能力。接觸剛度模型中，常用的是線性接觸模型，利用接觸的法向剛度Kn和切向剛度Ks兩個(gè)參數(shù)將顆粒間的法向、切向力和重疊量聯(lián)系起來(lái)。

(1)

(2)

其中：ξ={n,s}，n,s分別表示法向和切向?；颇Ｐ涂梢员硎緸?/p>

(3)

對(duì)于PFC2D中的接觸剛度模型和平行粘結(jié)模型，都可以視為兩顆粒及連接其圓心的彈性梁組成的系統(tǒng)[12]。設(shè)彈性梁的彈性模量為Ec(接觸彈性模量)，由彈性力學(xué)解可得

kn=2Ect。

(4)

對(duì)于混凝土材料而言，目前并沒有明確的標(biāo)準(zhǔn)規(guī)定采用哪些細(xì)觀本構(gòu)模型來(lái)進(jìn)行計(jì)算，表1列出的是不同學(xué)者研究混凝土問題時(shí)所采用的細(xì)觀模型參數(shù)。

表1 某些材料的細(xì)觀本構(gòu)模型參數(shù)取值Tab.1 Microparameters of micro-constitutive model for some materials

注:n_bond和s_bond分別為接觸粘結(jié)模型抗拉和抗剪強(qiáng)度

2 混凝土細(xì)觀本構(gòu)模型參數(shù)估計(jì)方法

巴西劈裂試驗(yàn)主要用來(lái)測(cè)試巖石材料的抗拉強(qiáng)度，通常是根據(jù)彈性力學(xué)解由試樣的最大承載力求得其抗拉強(qiáng)度。試樣在破壞前處于彈性變形階段，其荷載-位移曲線呈一近似直線。本研究主要目的是尋找一組合適的細(xì)觀參數(shù)使仿真得到的試樣宏觀力學(xué)響應(yīng)與試驗(yàn)結(jié)果吻合，因此重點(diǎn)關(guān)注試樣的最大承載力Fmax和其對(duì)應(yīng)的最大位移dmax。宏觀響應(yīng)與細(xì)觀參數(shù)之間存在著一定的函數(shù)關(guān)系，通常情況下，這些函數(shù)關(guān)系難以顯式表達(dá)出來(lái)，因此也就無(wú)法直接利用宏觀響應(yīng)確定細(xì)觀參數(shù)取值。為簡(jiǎn)化這些函數(shù)關(guān)系，假定宏觀響應(yīng)是細(xì)觀參數(shù)的線性疊加，且不同細(xì)觀參數(shù)之間不存在耦合，即：

(5)

(6)

其中：Xi表示細(xì)觀參數(shù)，ai和bi表示線性回歸系數(shù)，a0和b0為待定常數(shù)。為進(jìn)行細(xì)觀參數(shù)的估計(jì)，先給定一組細(xì)觀參數(shù)初始值，對(duì)每一個(gè)細(xì)觀參數(shù)，在其他參數(shù)取值不變的情況下，單獨(dú)改變?cè)搮?shù)的取值。分別取X={…,Xi,…}，X={…,Xi±ΔXi,…}和X={…,Xi±2ΔXi,…}5組細(xì)觀參數(shù)進(jìn)行PFC2D模擬，得到5個(gè)與Xi不同取值對(duì)應(yīng)的Fmax的值。以Xi為自變量對(duì)Fmax進(jìn)行線性回歸，得到回歸直線方程和相關(guān)系數(shù)R2，回歸直線的斜率即為Xi對(duì)Fmax的線性回歸系數(shù)ai。同理求出細(xì)觀參數(shù)對(duì)dmax的線性回歸系數(shù)bi。此外，通常認(rèn)為相關(guān)系數(shù)R2≥0.7時(shí)，兩變量具有相關(guān)性?；谶@一考慮，當(dāng)某一細(xì)觀參數(shù)對(duì)Fmax或dmax的相關(guān)系數(shù)R2<0.7時(shí)，認(rèn)為其對(duì)應(yīng)的線性回歸系數(shù)ai或bi為0，即這一細(xì)觀參數(shù)的取值變化對(duì)最大承載力Fmax或其對(duì)應(yīng)的最大位移dmax沒有影響。將ai和bi與對(duì)應(yīng)的細(xì)觀參數(shù)初始值相乘，所得結(jié)果與模擬得到的宏觀響應(yīng)之差即為待定常數(shù)a0和b0。

(7)

圖1 巴西劈裂試驗(yàn)PFC2D模型Fig.1 PFC2D model of Brazilian test

3 算例分析

翟越等[18]對(duì)C35商品混凝土試件了進(jìn)行巴西劈裂試驗(yàn)，試樣厚度50 mm，直徑100 mm。在PFC2D中建立這一試樣的離散元模型(如圖1)，試樣共由1 078個(gè)顆粒組成，顆粒半徑在1.0～2.0 mm內(nèi)均勻分布，顆粒厚度為50 mm，利用剛性墻對(duì)試樣施加荷載。根據(jù)文獻(xiàn)[14-15]，表2給出了細(xì)觀參數(shù)的初始值。

表2 C35混凝土細(xì)觀參數(shù)的初始值Tab.2 Initial estimation of C35 concrete microparameters

由圖2可以看出，隨著接觸彈性模量Ec的增大，試樣的最大承載力呈減小趨勢(shì)，最大承載力對(duì)應(yīng)的最大位移也逐漸減小。由圖3可以看出kn/ks取值的改變主要影響試樣最大承載力對(duì)應(yīng)的最大位移，最大位移隨kn/ks的增大而增大。kn/ks取值的改變對(duì)最大承載力的影響并無(wú)明顯規(guī)律。

圖2 Ec對(duì)荷載-位移曲線的影響Fig.2 Influence of Ec on load-displacement curves

圖3 kn/ks對(duì)荷載-位移曲線的影響Fig.3 Influence of kn/kson load-displacement curves

由圖4可以看出，隨著摩擦系數(shù)μ的增大，試樣的最大承載力也逐漸增大，最大承載力對(duì)應(yīng)的最大位移呈減小趨勢(shì)。摩擦系數(shù)的增大，限制了顆粒間滑移的發(fā)生，從而提高了承載力，減小了最大位移。

圖4 μ對(duì)荷載-位移曲線的影響Fig.4 Influence of μ on load-displacement curves

圖對(duì)荷載-位移曲線的影響Fig.5 Influence of on load-displacement curves

圖對(duì)荷載-位移曲線的影響Fig.6 Influence load-displacement curves

圖7 宏觀響應(yīng)隨Ec的變化Fig.7 Variation of macroscopic response versus Ec

根據(jù)圖2～6所獲得的信息，繪制不同細(xì)觀參數(shù)對(duì)應(yīng)的最大承載力和最大位移散點(diǎn)圖(圖7～11)。以各個(gè)細(xì)觀參數(shù)為自變量對(duì)最大承載力和最大位移進(jìn)行線性回歸，得到回歸直線方程和相關(guān)系數(shù)R2也顯示于圖7～11中。

圖8 宏觀響應(yīng)隨kn/ks的變化Fig.8 Variation of macroscopic response versuskn/ks

圖9 宏觀響應(yīng)隨μ的變化Fig.9 Variation of macroscopic response versus μ

圖10 宏觀響應(yīng)隨的變化Fig.10 Variation of macroscopic response

圖11 宏觀響應(yīng)隨的變化Fig.11 Variation of macroscopic response

表3中列出的是根據(jù)圖7～11中的回歸直線方程和相關(guān)系數(shù)確定的各細(xì)觀參數(shù)的線性回歸系數(shù)。將表2中C35混凝土細(xì)觀參數(shù)初始估計(jì)值及表3中相應(yīng)的線性回歸系數(shù)分別代入式(5)～(6)，即可計(jì)算得到待定常數(shù)ai和bi，從而得到宏觀響應(yīng)與細(xì)觀參數(shù)近似的定量關(guān)系(式(8)～(9))。

表3 細(xì)觀參數(shù)對(duì)宏觀響應(yīng)的線性回歸系數(shù)Tab.3 Linear regression coefficients of microparameters to macroscopic response

(8)

(9)

表4 C35混凝土細(xì)觀參數(shù)的估計(jì)值Tab.4 Estimated values of C35 concrete microparameters

圖12 PFC2D模擬與試驗(yàn)荷載-位移曲線對(duì)比Fig.12 Comparison of load-displacement curves between PFC2D simulation and experiment

圖12是使用細(xì)觀參數(shù)初始值和估計(jì)值進(jìn)行PFC2D模擬所得荷載-位移曲線與試驗(yàn)結(jié)果的對(duì)比，通過(guò)模擬得到的兩條曲線與試驗(yàn)曲線的相關(guān)系數(shù)分別為R2=0.767,R2=0.978,可以看出，使用估計(jì)值模擬得到的結(jié)果更加接近試驗(yàn)曲線。但這一曲線在接近峰值承載力時(shí)出現(xiàn)下降然后略有上升的現(xiàn)象，這是由于與壓板接觸處的顆粒因應(yīng)力集中而提前破壞后進(jìn)一步壓實(shí)所產(chǎn)生的。試樣斷裂面及PFC2D模擬裂紋分布情況分別如圖13(a)、(b)所示，圖13(b)中深黑色表示顆粒間膠結(jié)體拉伸破壞，灰色表示膠結(jié)體剪切破壞。

圖13 試驗(yàn)[18]與PFC2D模擬裂紋分布情況對(duì)比Fig.13 Comparison of crack distributions between test observation and PFC2D simulation

表5給出的是使用細(xì)觀參數(shù)初始值和估計(jì)值進(jìn)行模擬與試驗(yàn)曲線在不同位移處所得相對(duì)誤差的對(duì)比?？梢钥闯?，使用經(jīng)過(guò)優(yōu)化之后的細(xì)觀參數(shù)估計(jì)值計(jì)算得到的結(jié)果誤差明顯小于使用初始值得到的結(jié)果。

表5 使用不同參數(shù)在各個(gè)位移處的誤差對(duì)比Tab.5 Comparison of errors at each displacement using different parameters

由圖12～13和表5可以看出，使用表4中的細(xì)觀參數(shù)估計(jì)值，可以有效地模擬出混凝土試件在巴西劈裂試驗(yàn)中的力學(xué)行為和破壞過(guò)程，模擬曲線與試驗(yàn)曲線在整體變化趨勢(shì)上一致，最大承載力及其對(duì)應(yīng)的最大位移與試驗(yàn)結(jié)果比較接近。

4 結(jié)論

1) 離散元單元法的原理決定了只可以對(duì)顆粒的細(xì)觀參數(shù)進(jìn)行賦值，如何使細(xì)觀參數(shù)計(jì)算結(jié)果與宏觀試驗(yàn)結(jié)果吻合對(duì)PFC方法的應(yīng)用是極其重要的?；诨炷猎嚇拥陌臀髋言囼?yàn)，提出了使用一組試驗(yàn)數(shù)據(jù)并進(jìn)行少量PFC2D計(jì)算以得到一組滿足工程需求的細(xì)觀參數(shù)的方法，避免了試錯(cuò)法進(jìn)行參數(shù)調(diào)整的盲目性，節(jié)約了參數(shù)估計(jì)的時(shí)間。

2) PFC2D混凝土試樣的最大承載力與Ec和μ相關(guān)，最大承載力隨著Ec的增大而減小，隨著μ的增大而增大；試驗(yàn)最大承載力所對(duì)應(yīng)的最大位移與Ec和kn/ks相關(guān)，最大位移隨著Ec的增大而減小，隨著kn/ks的增大而增大；試驗(yàn)荷載-位移曲線的斜率也與Ec和μ相關(guān)，斜率隨著二者的增大而增大。

4) 利用本研究提出的細(xì)觀參數(shù)估計(jì)方法，得到了C35混凝土試件離散元細(xì)觀參數(shù)，PFC2D模擬結(jié)果表明，模擬曲線與試驗(yàn)曲線吻合較好，兩者相關(guān)系數(shù)R2為0.978，驗(yàn)證了該方法的有效性。

參考文獻(xiàn)：

[1]BIENIAWSKI Z T，HAWKES I.Suggested methods for determining tensile strength of rock materials[J].International Journal of Rock Mechanics & Mining Science & Geomechanics Abstracts,1978,15:99-103.

[2]ASTM C496-85.Standard test method for splitting tensile strength of cylindrical concrete specimens[S].Annual Book of ASTM Standards 4,1986:337-342.

[3]CHEN X D,GE L M,ZHOU J K,et al.Dynamic Brazilian test of concrete using split Hopkinson pressure bar[J].Materials and Structures,2017,50(1):1-15.

[4]CARMONA S,AGUADO A.New model for the indirect determination of the tensile stress-strain curve of concrete by means of the Brazilian test[J].Materials and Structures,2012,45(10):1473-1485.

[5]OHTSU M,KAWASAKI Y.AE-SiGMA analysis in Brazilian test and accelerated corrosion test of concrete[J].Journal of Acoustic Emission,2010,28:204-214.

[6]ABRISHAMBAF A,BARROS J,CUNHA V.Relation between fibre distribution and post-cracking behaviour in steel fibre reinforced self-compacting concrete panels[J].Cement & Concrete Research,2013,51(9):57-66.

[7]QIN C,ZHANG C H.Numerical study of dynamic behavior of concrete by meso-scale particle element modeling[J].International Journal of Impact Engineering,2011,38(12):1011-1021.

[8]GHAZVINIAN A,SARFARAZI V,SCHUBERT W,et al.A study of the failure mechanism of planar non-persistent open joints using PFC2D[J].Rock Mechanics and Rock Engineering,2012,45(5):677-693.

[9]MECHTCHERINE V,GRAM A,KRENZER K,et al.Simulation of fresh concrete flow using discrete element method (DEM):Theory and applications[J].Materials and Structures,2014,47(4):615-630.

[10]PAN J W,FENG Y T,JIN F,et al.Meso-scale particle modeling of concrete deterioration caused by alkali-aggregate reaction[J].International Journal for Numerical & Analytical Methods in Geomechanics,2013,37(16):2690-2705.

[11]KAZERANI T,YANG Z Y,ZHAO J.A discrete element model for predicting shear strength and degradation of rock joint by using compressive and tensile test data[J].Rock Mechanics and Rock Engineering,2012,45(5):695-709.

[12]ITASCA Consulting Group Incorporated.Particle flow in two-dimensions version 3.1 user manual[M].Minneapolis:ITASCA Consulting Group Incorporated,2005.

[13]張寧.碾壓混凝土碾壓特性的三維離散元數(shù)值模擬[D].杭州:浙江大學(xué),2015.

[14]張正珺,劉軍,胡文,等.混凝土材料破壞過(guò)程的二維離散元模擬[J].水力發(fā)電學(xué)報(bào),2010,29(5):22-27.

ZHANG Zhengjun,LIU Jun,HU Wen,et al.Two-dimensional simulation of concrete material fracturing by discrete element method[J].Journal of Hydroelectric Engineering,2010,29(5):22-27.

[15]唐陽(yáng),陶承品,徐霖玉,等.基于PFC2D的植生混凝土虛擬試驗(yàn)研究[J].水利科技與經(jīng)濟(jì),2012,20(12):10-13.

TANG Yang,TAO Chengpin,XU Linyu,et al.Research on virtual tests of green-growing concrete based on PFC2D[J].Water Conservancy Science and Technology and Economy,2012,20(12):10-13.

[16]叢宇.卸荷條件下演示破壞宏細(xì)觀機(jī)理與地下工程設(shè)計(jì)計(jì)算方法研究[D].青島:青島理工大學(xué),2014.

[17]李守巨,李德,于申.基于宏觀實(shí)驗(yàn)數(shù)據(jù)的堆石料細(xì)觀本構(gòu)模型參數(shù)反演[J].山東科技大學(xué)學(xué)報(bào)(自然科學(xué)版),2015,34(5):20-26.

LI Shouju,LI De,YU Shen.Meso-parameter inversion of constitutive model for rockfill materials based on marco experimental data[J].Journal of Shandong University of Science and Technology (Natural Science),2015,34(5):20-26.

[18]翟越,趙均海,艾曉芹,等.基于統(tǒng)一強(qiáng)度理論的巴西圓盤劈裂強(qiáng)度分析[J].建筑科學(xué)與工程學(xué)報(bào),2015,32(3):46-51.

ZHAI Yue,ZHAO Junhai,AI Xiaoqin,et al.Analysis on splitting of Brazilian disc based on unified strength theory[J].Journal of Architecture and Civil Engineering,2015,32(3):46-51.

[19]劉運(yùn)思,傅鶴林,饒軍應(yīng),等.不同層理方位影響下板巖各向異性巴西圓盤劈裂試驗(yàn)研究[J].巖石力學(xué)與工程學(xué)報(bào),2012,31(4):785-791.

LIU Yunsi,FU Helin,RAO Junying,et al.Research on Brazilian disc splitting tests for anisotropy of slate under influence of different bedding orientations[J].Chinese Journal of Rock Mechanics and Engineering,2012,31(4):785-791.