基于線性接觸模型的小麥三軸試驗(yàn)細(xì)觀模擬

曾長女，于 航

（河南工業(yè)大學(xué) 土木建筑學(xué)院，河南 鄭州 450001）

0 引言

在筒倉設(shè)計(jì)中，筒倉卸料引起的側(cè)壓力規(guī)律是關(guān)鍵要素之一.倉內(nèi)儲料的特性對筒倉側(cè)壓力具有較大的影響，其宏觀特性取決于儲料顆粒的細(xì)觀特性.利用數(shù)值試驗(yàn)，可方便地獲得需要的儲料特性計(jì)算參數(shù).離散元法用來處理散粒體流動特性具有很大的優(yōu)勢.PFC3D是ITASCA 公司基于Cundall 等[1]提出的離散元方法，是模擬不連續(xù)物體運(yùn)動及相互作用的三維顆粒離散元模擬程序.該程序被廣泛地應(yīng)用于模擬大變形、散體流動特性等研究中.

土工三軸試驗(yàn)的離散元模擬開展得較多.羅勇等[2]利用PFC3D對無黏性土的室內(nèi)常規(guī)三軸試驗(yàn)及其剪切帶形成和發(fā)展進(jìn)行了數(shù)值模擬.耿麗等[3]以PFC3D為工具，提出了考慮圍壓影響的粗粒土三軸試驗(yàn)細(xì)觀模型，并研究了不同的顆粒細(xì)觀參數(shù)對粗粒土強(qiáng)度的影響.邵磊等[4]通過生成不規(guī)則形狀的數(shù)值顆粒，模擬三軸排水剪切試驗(yàn).李識博等[5]利用PFC3D對隴西地區(qū)黃土三軸固結(jié)不排水剪切試驗(yàn)進(jìn)行模擬，考慮了不同圍壓、不同含水率對應(yīng)力、應(yīng)變等的影響.Cheng 等[6]運(yùn)用PFC3D對巖石的團(tuán)塊組合進(jìn)行模擬，研究了三軸應(yīng)力-應(yīng)變關(guān)系.Hong-hua Zhao 等[7]采用PFC3D中的線性接觸模型，研究了摩擦系數(shù)對三軸試驗(yàn)獲得的密砂應(yīng)力-應(yīng)變的影響.Backstrom 等[8]利用顆粒離散元研究了花崗巖細(xì)觀參數(shù)及其應(yīng)力-應(yīng)變曲線.現(xiàn)有的研究多集中于不同土體三軸試驗(yàn)，基于糧食散體細(xì)觀參數(shù)的研究很少，而利用離散元法分析筒倉卸料時(shí)，細(xì)觀參數(shù)多直接設(shè)定，沒有更深入地說明.因此有必要獲得小麥等糧食散體的細(xì)觀參數(shù)，為今后系統(tǒng)地研究糧食散體的細(xì)觀行為提供基礎(chǔ)參數(shù).

小麥的宏觀特性必然取決于其顆粒的細(xì)觀特性，離散元法從細(xì)觀角度對小麥的力學(xué)特性進(jìn)行數(shù)值模擬，分析小麥的細(xì)觀參數(shù)變化對小麥的宏觀力學(xué)行為的影響，再現(xiàn)小麥剪切試驗(yàn)曲線.在進(jìn)行數(shù)值計(jì)算和理論研究時(shí)，小麥堆體的彈性模量等參數(shù)的確定對筒倉結(jié)構(gòu)設(shè)計(jì)具有重要意義.作者基于小麥三軸試驗(yàn)，通過PFC3D模擬室內(nèi)小麥三軸壓縮試驗(yàn)，經(jīng)反復(fù)試算，標(biāo)定細(xì)觀參數(shù)，并反算獲得其宏觀參數(shù).

1 PFC3D 模擬過程

1.1 細(xì)觀計(jì)算模型的選取

PFC3D中的接觸模型包括3 個(gè)模型：接觸剛度模型；滑移-分離模型；黏結(jié)模型.接觸剛度模型描述了接觸力和相對位移的彈性關(guān)系如圖1 所示，采用線性接觸剛度模型，該模型應(yīng)用較多，細(xì)觀參數(shù)較少，也可滿足模擬結(jié)果的精度要求.

在線性接觸剛度模型中，接觸剛度為不變值，應(yīng)力和相對位移是線性相關(guān)的.法向作用力滿足式（1）：

式中：Kn為法向剛度，屬割線模量；Un為法向位移.

切向作用力滿足式（2）：

式中：Ks為切向剛度；ΔUs為切向位移增量.

圖1 線性接觸剛度模型Fig.1 Linear contact stistness model

在細(xì)觀模擬中，有顆粒法向剛度Kn、顆粒切向剛度Ks、顆粒間摩擦系數(shù)μ、三軸試驗(yàn)加載速度v等細(xì)觀參數(shù)需重點(diǎn)考慮.

1.2 三軸細(xì)觀模型的建立

1.2.1 生成孔隙率為n 的試樣

通過顆粒半徑放大法使模型試樣達(dá)到目標(biāo)孔隙率.對所有顆粒采用相同的半徑放大系數(shù)m，即：

式中：n0為初始孔隙率；n 為目標(biāo)孔隙率；R0為初始顆粒半徑.

模型試樣生成時(shí)，先假定半徑放大系數(shù)m，通過式（3）反算出初始孔隙率：

然后設(shè)定顆粒半徑分布范圍是（Rmin，Rmax），顆粒平均半徑為：

試樣尺寸為h×d，則生成的顆粒數(shù)目為：

根據(jù)前面步驟反復(fù)計(jì)算，直至達(dá)到設(shè)定的孔隙率.

1.2.2 墻體剛度的選取

生成的數(shù)值模型如圖2 所示，包括圓柱形墻體和上、下兩面墻體.圓柱形墻體模擬三軸試驗(yàn)中的橡皮膜，上、下兩面墻體模擬三軸試驗(yàn)中的加載墻體.墻體剛度的選擇對試驗(yàn)結(jié)果影響較大.已有研究表明，墻體剛度過大會導(dǎo)致不容易平衡，太小則導(dǎo)致側(cè)壓力過小，從而與實(shí)際不符.在數(shù)值模擬中，可采用圓柱形墻體的剛度為顆粒剛度0.1 倍時(shí)，可很好地模擬橡皮膜的柔性約束，上、下墻體的剛度為顆粒剛度10 倍以上時(shí)，模擬剛性加載面效果更好.因此，作者選擇圓柱形墻體的剛度為顆粒剛度的0.1 倍，上、下加載墻體的剛度為顆粒剛度的10 倍[9].

圖2 三軸試驗(yàn)試樣數(shù)值模型Fig.2 Numerical model of triaxial test sample

1.2.3 三軸試驗(yàn)加載模擬

三軸試驗(yàn)?zāi)Ｐ图虞d時(shí)，圓柱形墻體用來模擬橡皮膜，可施加力模擬三軸試驗(yàn)中的圍壓.通過在上、下兩墻施加力模擬三軸試驗(yàn)的軸向壓力.具體計(jì)算公式如式（8）、（9）、（10）.

式中：r 為試樣底面半徑.

式中：h 為試樣高度.

軸向應(yīng)變和徑向應(yīng)變都可用下式計(jì)算：

式中：L 為當(dāng)前式樣的高度或?qū)挾?；L0為初始試樣的高度或?qū)挾?

試驗(yàn)加載步驟為：

（1）通過伺服機(jī)制給圓柱墻一個(gè)速度控制圍壓不變.

（2）給上、下墻大小相同、方向相反的速度壓縮試樣：

式（13）、（14），以使模型達(dá)到目標(biāo)的圍壓和軸向荷載.

（4）計(jì)算參數(shù)G 的選擇.

在一個(gè)計(jì)算時(shí)步中，墻體應(yīng)力的增量為：

式中：Nc為 顆粒與墻體的接數(shù)目；為所有和墻體接觸的平均剛度；A 為墻體的面積.

求解過程需要一定的穩(wěn)定性，為確保加載過程的穩(wěn)定，在實(shí)際操作中增加一個(gè)松弛因子α（一般為0.5）.

由式（11）、（12）、（15）和（16）可得：

為了獲得穩(wěn)定條件，把參數(shù)G 定義為：

在每一次循環(huán)開始時(shí)，墻的速度給定一個(gè)值滿足式（11）、（12），而G 要滿足式（18）.

2 細(xì)觀模擬結(jié)果

PFC3D三軸模型參數(shù)的篩選與輸入，是一個(gè)非常費(fèi)時(shí)的過程，生成的顆粒數(shù)量對計(jì)算速度影響很大，經(jīng)過多次試算，本次試驗(yàn)?zāi)M采用的試樣尺寸取值0.30 m×0.15 m，顆粒半徑取值范圍定為（0.004，0.008）m，孔隙率為0.3，該模型生成3 877個(gè)顆粒.

2.1 應(yīng)力-應(yīng)變曲線

根據(jù)已獲得的室內(nèi)試驗(yàn)的應(yīng)力-應(yīng)變曲線，不斷調(diào)整細(xì)觀參數(shù)并反復(fù)試算，模擬得出的應(yīng)力-應(yīng)變曲線與真實(shí)試驗(yàn)得出的應(yīng)力-應(yīng)變曲線對比如圖3 所示.由圖3 可知，試驗(yàn)曲線與數(shù)值模擬曲線擬合較好，由此可獲得合適的細(xì)觀參數(shù)，具體如表1 所示.

圖3 試驗(yàn)與數(shù)值模擬的應(yīng)力-應(yīng)變曲線Fig.3 Stress-strain curve from test and numerital simulation

表1 材料參數(shù)Table 1 Material parameters used in simulations

2.2 三軸試驗(yàn)中小麥宏觀參數(shù)推算

（1）彈性模量E

E 的確定如圖4 所示，可計(jì)算獲得小麥的宏觀參數(shù).對于三軸試驗(yàn)而言，采用式（19）可獲得其彈性模量.初始階段應(yīng)力-應(yīng)變曲線的斜率可取為彈性模量：

圖4 E 的確定Fig.4 Stress-strain curve from simulation

（2）泊松比ν

由試驗(yàn)確定泊松比比較困難，通過如圖5 所示體變-軸向應(yīng)變曲線可確定泊松比.

（3）內(nèi)摩擦角φ

該模擬試驗(yàn)時(shí)圍壓σ3保持不變，偏應(yīng)力最大值σ1max，由圖4 可知，利用公式（20）可計(jì)算得到顆粒的內(nèi)摩擦角：

圖5 泊松比確定Fig.5 Poisson ratio

（4）剪脹角ψ

體變-軸向應(yīng)變曲線與x 軸有一交點(diǎn)，利用最小二乘法擬合出與該交點(diǎn)往后的曲線最接近的直線，由公式（21）可求出剪脹角ψ：

2.3 與文獻(xiàn)結(jié)果對比

利用前面計(jì)算方法，小麥的宏觀參數(shù)計(jì)算結(jié)果如表2 所示.為了驗(yàn)證模擬結(jié)果，對比了文獻(xiàn)中關(guān)于小麥的彈性模量等宏觀參數(shù).文中試驗(yàn)?zāi)M的是100 kPa 圍壓下的應(yīng)力應(yīng)變關(guān)系，由表2 可知，小麥的宏觀參數(shù)如彈性模量等值差異性較大，一方面，采用不同測試方法測定會導(dǎo)致數(shù)值誤差；另一方面，小麥本身的硬度、含水量等初始狀態(tài)都會對彈性模型產(chǎn)生較大的影響.文中獲得的彈性模量值在文獻(xiàn)的數(shù)值范圍之內(nèi).

表2 小麥的宏觀參數(shù)對比Table 2 Comparison of wheat maoropameters

3 結(jié)論

利用PFC3D離散元方法模擬室內(nèi)小麥三軸試驗(yàn)，獲得小麥三軸試驗(yàn)的應(yīng)力-應(yīng)變曲線，通過對比宏觀試驗(yàn)曲線，標(biāo)定了小麥的細(xì)觀參數(shù)，得出以下結(jié)論：

（1）PFC3D數(shù)值模擬能靈活改變?nèi)S試驗(yàn)的邊界、初始條件等，可有效地模擬室內(nèi)試驗(yàn)，是進(jìn)行三軸試驗(yàn)的有效數(shù)值手段.

（2）數(shù)值試驗(yàn)結(jié)果表明，小麥顆粒材料的宏觀參數(shù)彈性模量、泊松比、內(nèi)摩擦角和剪脹角與細(xì)觀參數(shù)顆粒剛度、顆粒間摩擦系數(shù)等具有良好的相關(guān)性.

（3）通過數(shù)值試驗(yàn)，獲得小麥的細(xì)觀參數(shù)，對比分析得到的宏觀參數(shù)符合要求.

［1］Cundall P A，Strack O D L.A discrete numerical model for granular assemblies［J］.Geote -chnique，1979，29（1）:47-65.

［2］羅勇，龔曉南，連峰.三維離散顆粒單元模擬無黏性土的工程力學(xué)性質(zhì)［J］.巖土工程學(xué)報(bào)，2008，30（2）：292-297.

［3］耿麗，黃志強(qiáng)，苗雨.粗粒土三軸試驗(yàn)的細(xì)觀模擬［J］.土木工程與管理學(xué)報(bào)，2011，28（4）：24-29.

［4］邵磊，遲世春，張勇，等.基于顆粒流的堆石料三軸剪切試驗(yàn)研究［J］.巖土力學(xué)，2013，34（3）:711-720.

［5］李識博，王常明，王念秦，等.黃土三軸試驗(yàn)的顆粒流數(shù)值模擬［J］.中國公路學(xué)報(bào)，2013，26（6）：22-29.

［6］Cheng Y P，Nakata Y，Bolton M D.Discrete element simulation of crushable soil ［J］.Geotechnique，2003，53（7）：633-641.

［7］Zhao Hong-hua，Shao Long-tan，Ji Shunying.Numerical simulation of triaxial test on the dense sand by DEM［J］.Geotechnical Special Publication，2011，220：233-241.

［8］Backstrom A，Antikainen J，Backers T.Numerical modelling of uniaxial compressive failure of granite with and without saline porewater［J］.International Journal of Rock Mechanics&Mining Sciences，2008，45（7）:1126-1142.

［9］羅勇.土工問題的顆粒流數(shù)值模擬及應(yīng)用研究［D］.杭州：浙江大學(xué)，2007.

［10］程緒鐸，安蓉蓉，曹陽，等.小麥糧堆彈性模量的實(shí)驗(yàn)測定與研究［J］.糧食儲藏，2009，38（6）：22-25.

［11］程緒鐸，陸琳琳，石翠霞.小麥摩擦特性的試驗(yàn)研究［J］.中國糧油學(xué)報(bào)，2012，27（4）：15-19.

［12］石翠霞.筒倉內(nèi)小麥、玉米堆的壓縮特性、倉壁壓應(yīng)力及儲量總重量的研究［D］.南京：南京財(cái)經(jīng)大學(xué)，2012.

［13］梁醒培，楊伯源，王西院.深筒倉動態(tài)壓力數(shù)值模擬［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2008，31（3）：427-429.