土工袋双轴压缩的离散元数值模拟

2019-03-11高军军刘斯宏王恩准

水利水电科技进展 2019年1期

高军军,汤雷,刘斯宏,李建,王恩准

(1.南京水利科学研究院材料结构研究所,江苏南京 210029;2.河海大学水利水电学院,江苏南京 210098; 3.南京瑞迪高新技术有限公司,江苏南京 210029)

土工袋技术是近年来通过系列试验研究、理论分析及工程应用而发展起来的一项地基处理新技术。Matsuoka等[1-2]通过一系列的室内试验研究并验证了土工袋的加筋机理,即将散状土体装入土工编织袋,利用压缩过程中袋子产生的张力约束袋内土体,达到提高袋内土体强度的目的。土工袋已在房屋基础、公路与铁路路基、沟槽回填及挡土墙等工程中得到广泛应用。目前对土工袋的研究主要以室内试验为主：刘斯宏等[3]通过现场平板载荷试验验证了土工袋具有比一般土体高得多的承载能力;高军军等[4-5]通过室内试验和数值模拟验证土工袋可显著提高淤泥土的强度,其极限抗压强度约为普通混凝土的1/10;李玲君等[6]探讨了不同袋内充填材料对土工袋动力特性的影响;相关学者还通过铝棒摩擦试验研究了不同条件下柔性土工袋层间摩擦与刚性体间摩擦特征的区别[7]。

对土工袋的数值模拟研究相对较少,刘斯宏等[8]将袋子张力的作用等效为附加应力,利用弹塑性有限元法对加固后的软基承载力进行了数值模拟;Ansari等[9]将土和袋子分别采用不同的模型利用有限元法研究单个土工袋在竖向荷载下的力学反应。但土工袋是一种复杂的复合型材料,存在非均质性和介质不连续性,尤其是不同材料介质间的接触问题,采用基于连续介质理论的有限元法反映其真实力学特性的难度较大,土工袋固有的间隙也直接制约有限元法计算的有效进行。故而采用基于非连续介质理论的颗粒离散元法对土工袋的力学特性进行数值模拟来验证土工袋的加筋机理有一定的现实意义。本文对土工袋在双轴压缩条件下的力学特性进行离散元数值模拟,从细观角度分析土工袋内部应力分布及发展规律。

1 土工袋模型的生成

考虑真实土体颗粒分布的空间随机性,本文采用蒙特卡洛随机算法生成一定颗粒粒径及颗粒级配的袋内土体样本,其中土颗粒间采用弹簧-阻尼器-滑块接触模型,袋子颗粒间采用皮筋-阻尼器接触模型。

1.1 蒙特卡洛随机算法

蒙特卡洛随机算法又称计算机随机模拟方法,它是一种基于随机数的计算方法[10-12],在大数定理和中心极限定理等概率论的基础上,对随机变量进行方差判定来分析随机问题的不确定程度。

设X1、X2、…、Xn(n为随机变量的总数)为独立且均匀分布的随机变量序列,其期望值为E,标准差为S,则正态随机数可表示为

(1)

颗粒生成程序中,根据边界条件和颗粒数可以求出颗粒的平均半径,输入比平均半径较小的一个值R0和预设标准差S0,通过式(2)就可以得到第i个颗粒粒径Ri,如果Ri小于0,则回到随机算法重新生成100个处在区间[0,1]内的随机变量并重新得到Q100,根据式(2)得到新的粒径值,重复判断,直到Ri大于0,进入下一个循环生成第i+1个颗粒的粒径Ri+1。

Ri=R0+S0Q100

(2)

根据颗粒单元的交叉(重合)情况通过分级衰减粒径大小的方法再重新定位颗粒单元,如此往复直至生成较满意的颗粒系统。此颗粒系统需满足两个条件:①颗粒在给定空间内具有各向均匀性;②颗粒体的初始重叠量不宜过大。

1.2 颗粒体生成

首先确定矩形试样尺寸为长度40 cm、高度10 cm(土工袋的合理尺寸[13])。如果颗粒最大粒径为9.0 mm,生成初始颗粒数为1 250,再根据衰减系数对发生重叠的颗粒进行分级粒径衰减,经过5级衰减后,得到粒径分别为9.00 mm、7.20 mm、6.12 mm、5.20 mm、4.16 mm、3.33 mm。此时生成的颗粒样本虽然不存在重叠现象,但是部分颗粒间隙较大,存在颗粒腾空现象,需要再进行自重堆积得到新的颗粒位置状态。此时颗粒样本在顶部并不平整,会影响竖向均布荷载的施加,可在其顶部先施加10 kPa的预荷载,同时限制侧向变形为0,待稳定后得到试验最终样本。

最终生成的样本包含1 250个颗粒,其中颗粒数按粒径从大到小分别为:180、148、156、189、184、393,颗粒级配曲线如图1所示,最终的颗粒孔隙比为0.20,不均匀系数为2.1。

图1 颗粒级配曲线

1.3 颗粒接触模型

在颗粒离散单元中,宏观系统的力学行为通过细观颗粒间的接触模型体现出来。基于牛顿第二定律,通过颗粒间的作用力求得各颗粒的加速度,加速度对时间进行一次积分即可求得颗粒的速度,加速度对时间进行二次积分则可求得颗粒的位移。

在土工袋的离散元模拟计算中,袋内土颗粒间的接触简化为弹簧-阻尼器-滑块模型,如图2(a)所示,弹簧主要模拟颗粒间的压缩和切向所用,滑块主要模拟颗粒间的滑移,而阻尼器主要是为了消耗系统中多余的能量使系统快速趋向平衡状态;土工袋通过在土颗粒的外面增加一层具有张力的小颗粒模拟,袋子颗粒间的接触简化为皮筋-阻尼器模型,如图2(b)所示,其间的接触模型与土颗粒相同,但只有法向接触,无切向接触,且只受拉不受压。

图2 袋子DEM数值模拟颗粒接触模型

2 无袋子情况下的土体强度

为了反映土工袋对袋内土体的加筋作用(即提高袋内土体强度),先对无袋子情况下近似同等体积的土体进行数值模拟,分析其在双轴压缩条件下的力学强度。

表1 土体试样DEM计算参数

模拟试验在4块相对独立、互不影响的刚性加载板的范围内进行,如图3所示。颗粒密度为2.7 g/cm3,颗粒间的摩擦角为20°;假定颗粒与刚性加载板之间无摩擦,避免加载板对试样端部变形及试样内部应力的影响。具体计算参数如表1所示。

图3 无袋子情况下的土体试样双轴压缩模拟

加载过程中对试样施加以50 kPa的稳定围压,再分别对顶部和底部的加载板施加以大小为0.01 cm/s、方向相反的竖向等应变率加载速度。

图4为无袋子包围的土颗粒试样偏应力及体应变与轴向应变之间的关系曲线。从图中可以看出土体在双轴压缩过程中应力随着应变非线性增长。而土体的体应变在压缩过程中先发生微小的剪缩现象,随即开始剪胀,且体应变随轴向应变的增大而增大,当土体发生破坏后,剪胀值有所降低。

图4 无袋子土体的应力-应变关系曲线

图5 土颗粒试样内部力链分布

图5为土颗粒试样达到极限强度时的力链结构图,线条的粗细代表颗粒间接触力的大小。力链是颗粒物质力学的基本概念,孙其诚等[14]认为颗粒在外荷载作用下相互挤压形成接触网络,该网络是支撑外荷载的物质基础,就如同土力学中土骨架的概念一样[15-16]。从图5可看出,颗粒间的强力链(即较粗力链)分布较均匀,基本上与大主应力方向平行,而两侧颗粒间的力链较弱,说明土颗粒试样的破坏由竖向荷载主导,这与事实相符。

为了通过摩尔圆来分析土体的强度特性,又对该试样作了不同围压下(100 kPa、200 kPa)的双轴压缩模拟。从试验样本的摩尔强度线可以看出,不同围压条件下的应力摩尔圆基本相切于一直线,故本次试验样本的内摩擦角可取24.3°。

3 土工袋的力学强度

为了研究土工袋的力学性能,只需在原试验样本上加一圈袋子小颗粒,如图6所示,假设颗粒级配及样本内摩擦角不变。袋子颗粒间保留一定距离,允许袋子在外荷载作用下发生伸缩变形;且假设间距相等,为袋子颗粒粒径的2.5倍。该土工袋试样模型中,土颗粒数为1 199,袋子颗粒数为396,袋子颗粒粒径为0.1 cm,其密度设为0.8 g/m3。土颗粒的计算参数如表1所示,袋子颗粒的计算参数见表2,模拟中不考虑袋子与刚性边间的摩擦作用。

加载方式与无袋子情况相同。需要注意的是,在计算主应变时,将土工袋简化为理想的矩形结构,即利用上下刚性边的相对位移计算土工袋的轴向应变,利用左右刚性边的相对位移计算土工袋的侧向应变。

表2 土工袋试样DEM计算参数

图6 土工袋试样双轴压缩模拟

图7为土工袋在围压50 kPa条件下的偏应力-轴向应变和体应变-轴向应变关系曲线。从图7可看出,土工袋的偏应力随着轴向应变的增长逐渐增大,直到土工袋发生破坏后,迅速下降;而土工袋的体应变在加载初期也发生微小的剪缩现象,随后快速膨胀,到后期增长速率变小,这与无袋子情况下的土体变形规律一致。不同的是,土工袋破坏时的偏应力(680 kPa)远大于无袋子的情况(70 kPa),并且土工袋破坏时的总体变形明显大于无袋子的情况。袋子具有可延性,到后期袋子的延伸速率下降,逐渐达到袋子的抗拉强度。从图中还可以看出,随着轴向应变的增长,偏应力有好几处(红圈标记)发生一定的下降后重新上升,这是由于虽然土工袋仍未破坏,但此时袋内的土体早已达到极限状态,袋内土体产生局部应力集中,部分土颗粒发生旋转,在袋子的侧向约束下重新达到新的平衡状态,而在此过程中体应变并无明显变化。

图7 土工袋的应力-应变关系曲线

图8(a)为土工袋破坏前某一状态下(轴向应变ε1=3%时)相邻袋子颗粒间的张力大小沿土工袋周边的分布情况,可见,袋子张力分布并不是均匀的,在袋子不同的位置张力大小不同,但基本上在一个平均值附近上下波动。对所有袋子颗粒间的张力取平均值可以得到袋子平均张力的变化情况,图8(b)为模拟过程中袋子的平均张力随着轴向应变的变化。袋子的平均张力并不呈线性增长趋势,这是由于在变形过程中袋子周长的伸长与土工袋轴向应变并不呈线性关系。袋子破坏时袋子平均张力接近12.5 kN/m。

图8 土工袋袋子张力的变化

根据刘斯宏等[2]对土工袋加筋原理的阐述,将袋子张力等效为作用于土工袋轴向及侧向的附加应力(2T/B,2T/H),强度公式可表示为

(3)

图9 袋内土体主应力之间的关系

式中:σ1为土工袋大主应力;σ3为土工袋小主应力;T为袋子的张力;B为土工袋的宽度;H为土工袋的高度;KP为被动土压力系数,KP=(1+sinφ)/(1-sinφ),φ=24.3°。

取图7中红圈内的大小主应力,再结合图8中对应轴向应变下的平均张力,可以得到袋内土体主应力之间的关系如图9所示,大小主应力之间基本呈线性关系,且斜率近似与KP一致,验证了理论强度式(3)的合理性。

假设土工袋在压缩过程中体积不变,可以推得土工袋的轴向应变ε1与侧向应变ε3之间的关系为ε3=ε1/(ε1-1)。图10显示DEM数值模拟得到的土工袋轴向应变与侧向应变之间的关系(黑色实线),发现在加载初期土工袋基本上满足体积不变的假设,但当轴向应变超过2%以后,真实的侧向应变绝对值大于假设体积不变得到的值,即土工袋的体积发生了变化。DEM数值模拟中的ε1与ε3的关系更接近于多项式表达式(如图10红色虚线所示):

ε3= 0.007 85ε13-0.138 34ε12-

0.800 65ε1-0.049 87

(4)

图12 土工袋内部颗粒的细观分析

图10 土工袋轴向应变与侧向应变的关系

通过几何换算,袋子的应变εbag可以表示为

(5)

式中:m为初始状态下土工袋的宽度与高度的比值。

将式(5)代入T=kεbag(其中k为袋子接触模型的抗拉系数),并与式(4)一起代入式(3)可以推导出适合此双轴压缩试验的应力-应变关系:

(6)

式中:H0、B0分别为土工袋的初始高度和宽度。图11为DEM数值模拟结果与理论计算值的对比,可见推导公式与模拟结果吻合较好。

图11 DEM数值模拟结果与推导公式的对比

4 土工袋细观结构变化

图12为不同加载时刻的土工袋细观分析图。ε1=0时为围压50 kPa的情况,因为颗粒处于等向压缩状态,土工袋的内部应力比较均匀,力链群形成环状分布,表明此时外荷载由土颗粒均匀承担。当上下刚边以一定速度压缩土工袋时,土工袋被压扁,内部应力重新分布,ε1=2.8%时可以观察到内部应力出现局部集中化,力链方向大致为大主应力的方向;而土颗粒的位移较小且较乱,由于此时偏应力不大,在双轴压缩的情况下,土颗粒大致都往内部移动,导致土工袋在加载初期出现微小剪缩现象,同时袋子4个边角上的土颗粒有向外挤出的现象。当ε1=5.1%时,内部应力向中部集中,颗粒位移加大,颗粒基本上向土工袋两侧跑动,中部的颗粒位移较小;当ε1=11.1%时,土体内部出现滑移带,颗粒有明显向两侧发生大位移的趋势;随着轴向应力的继续增加,内部颗粒应力重新调整,直至土工袋内部基本上都处于强力链的状态,在轴向应变达到12.6%时发生破坏。在加载后期,颗粒的运动也呈对称性地向两端扩散,出现明显的滑移带,直至土工袋发生破坏时,土颗粒有明显向右下侧滑移的趋势(图12(i)中圆圈处),这是由于袋子就在此处发生了破坏(如图12(j)所示),对袋内土体失去约束作用,导致颗粒向外跑动,同时发现此时强力链开始从右下侧消失,表明袋子破坏是土工袋达到极限强度的直接原因。