1、第 39卷 第 1期 三峡大学学报（自然科学版） Vol.39 No.12017年 .1月 Journal of China Three Gorges University (Natural Sciences) Jan. 2017三维离散元大尺度模拟中变粒径方法作者 11，作者 22，1 ，作者 3*1，作者 41(1. 上海海事大学海洋科学与工程学院，上海 201306；2. 中交上海航道勘察设计研究院有限公司，上海 200120)摘 要：在已有的变粒径建模方法基础上，针对不同粒径区宏观弹性不一致和界面区大小颗粒互相渗入两个问题进行优化。通过在近场核心区域采用小粒径颗粒模拟真实土颗粒，远场边
2、界区域逐级放大颗粒粒径，不同粒径之间采用相同密实度的混合粒径模拟，实现三维大尺度离散元模拟，以减少模拟颗粒数量提升运算速度。最后，通过模拟静力触探试验（CPT），表明该方法的有效性。该方法可为三维离散元模拟工程大尺度问题提供新途径。关键词：离散元模拟；变粒径；弹性常数；尺度关系；相对密实度；静力触探中图分类号：TU43 文献标识码： A 文章编号：10004548(2017)01006209作者简介：作者 1(1980 )，男，讲师，博士，主要从事土体细观力学与离散元模拟等方面的教学和科研。 E-mail: Refinement and application of variable part
3、icle-size methods in 3D discrete element modelling for large-scale problemsAuthor11, Author22, 1, Author31, Author41(1. College of Ocean and Science Engineering, Shanghai Maritime University, Shanghai 201306, China; 2. Shanghai Waterway Engineering Design and Consulting Co., Ltd., Shanghai 200120, C
4、hina)Abstract: A refinement method is proposed to solve the difference in the macroscopic elasticity and the interpenetration in interface area between different particle-size areas in the DEM. In simulations of 3D DEM large-scale problems, small balls are used in the near field core area to simulat
5、e the real soil, while the far field boundary regions are filled with larger balls. In addition, a mixed area is set between different particle-size areas with the same relative density. The number of particles in the model can be reduced enormously, which improves largely the calculation efficiency
6、. For scaling particle systems, the scale relationship is established on the static condition by the identical macro elasticity of uniform granules. This study may provide a new approach to the application of 3D discrete element modelling in the large-scale geotechnical engineering.Key words: discre
7、te element modeling; variable particle-size; elastic constant; scaling law; relative density; CPT0 引 言砂土具有天然散粒特性，虽然传统的连续介质力学方法利用本构方程可以模拟一些工程问题，但对于局部的大变形破坏区域和土与结构非连续接触作用区域，连续性假设不再有效。离散元利用颗粒间的接触作用可从颗粒尺度上模拟散粒材料的力学行为，目前，成为研究岩土材料大变形破坏等问题的重要方法。早期学者们多采用二维圆盘离散元模拟 1-3，但二维圆盘运动约束不同于真实的颗粒。而三维离散元从颗粒尺度上模拟工程尺度问题往往
8、需要大量颗粒数目，如何提高计算效率成为学者们所关注的问题。1 变粒径建模方法在 DEM 模拟中，为了减小颗粒数目，将所模拟问题划分若干区域，如图 1 所示，从近场到远场划分为 Z1，Z 2， .，Z n区域，在近场核心区域采用小颗粒模拟，以模拟其复杂的物理力学特性；其它各区域模拟颗粒粒径逐级放大，以模拟边界。同时在相邻变粒径区域增设混合区域 M12，M 23， .，M (n-1)n，防止相邻粒径交界处出现相互嵌入混合。混合区为相邻区域的两种颗粒粒径混合，密实度与相邻区域一致。由于各区域粒径进行了放大，需要按照尺度关系调整颗粒材料参数，以保证不同区域相应的物理力基金项目：国家自然科学基金青年基金
9、项目（） ； 收稿日期： *通讯作者（E-mail: ）DOI：2 三峡大学学报（自然科学版） 2017 年2学特性的一致性。2 不同粒径颗粒体系的弹性验证2.1 粒径不同、细观结构相同试样生成方法将某一区域小颗粒体系放大至另一区域中，颗粒粒径放大 倍，按下式构建大颗粒体系坐标系：， (1)xc式中， 为坐标平移矢量。c在大颗粒区域生成与小颗粒相对应的大颗粒，如图 3 所示，粒径为小颗粒的 倍，位置矢量为相应小颗粒位置矢量的 倍，即(2)nd，x式中， 和 分别为大和小颗粒体系中第 n 个球的nx位置矢量。图 1 颗粒粒径放大示意图对于任意相接触的 a，b 两球，小、大颗粒中心连线矢量分别为(
10、3)baba ()c c， 。lxxl说明此时大小颗粒体系的接触分布一致，即接触组构相同。2.2 三轴试验模拟验证通过 DEM 模拟三轴试验，比较分析不同粒径颗粒体系的宏观弹性特性，验证不同粒径颗粒体系宏观弹性一致情况下的颗粒参数尺度关系。模拟的工况列于表 1 中。随机生成 3 种不同粒径（10，30，60 mm） 、孔隙比和尺寸相同的三轴试样（见图 4） 。#1， #2 和 #3 试验中，3 种试样采用相同的细观参数。#4 和 #5 试验中，保持刚度比 相同，按照法向刚度与粒径比与 #1 中相同，相应增大刚度。 #6 和 #7 试样按照上述粒径不同、接触组构相同试样生成方法，在#1 试样基础
11、上等比例放大，生成与 #1 组构相同试样，按照法向刚度与粒径比与 #1 中相同，相应增大刚度。施加围压 5 MPa 进行小应变弹性模拟，采用线弹性接触模拟接触作用，模拟结果见图 5。表 1 模拟工况编号粒径/mm法向刚度kn/(Nm-1)切向刚度ks/(Nm-1)摩擦系数 组构#1 10 8107 8105 0.8 #2 30 8107 8105 0.8 #3 60 8107 8105 0.8 #4 30 2.4108 2.4106 0.8 同 #2#5 60 4.8108 4.8106 0.8 同 #3#6 30 2.4108 2.4106 0.8 同 #1#7 60 4.8108 4.81
12、06 0.8 同 #1注： #1， #2， #3 为随机生成试样； #6， #7 与 #1 试样接触组构一致，法向刚度与粒径比保持相同。3 静力触探变粒径 DEM 模拟采用上述变粒径方法模拟静力触探，与全部小颗粒模拟结果对比分析，以验证变粒径方法的有效性。根据魏杰 20理论推导研究，塑性区域在 4D（D为探杆直径）范围内。Huang 等 21的离散元静力触探模拟研究表明探头对土体的应力应变响应主要也集中在 4D 范围以内。因此在本文的模型中小颗粒集合体区域半径为 75 mm4D（72 mm），模拟静力触探大变形，75150 mm 区域采用变粒径模拟。模拟分区情况和各区域采用的颗粒粒径列于表 3
13、，A 组、B组和 C 组为 McDowell 等 8结果， A 组采用均一小颗粒模拟，将其结果作为参照，B 组和 C 组分区分别采用 2 种和 3 种粒径模拟，但大小颗粒细观参数相同，法向接触刚度 kn=5105 N/m，切向接触刚度ks=5105 N/m，摩擦系数为 0.5，试样孔隙率为0.3545。D 组和 E 组采用 3.1 节所述方法生成组构与小颗粒区相同的大颗粒体系，同时，通过保持刚度比和法向刚度与粒径比相同，根据颗粒粒径，优化大颗粒刚度参数。D 组采用 3 种粒径模拟细观参数分别为kn=5105 N/m，k s=5105 N/m；k n=7.5105 N/m，k s=7.5105
14、N/m 和 kn=1106 N/m，k s=1106 N/m，E 组按 4.1 节所述方法在大小颗粒区中间设置一层粒径混合区（见图 13） ，其颗粒细观刚度为kn=5105 N/m，k s=5105 N/m 和 kn=1.5106 N/m，k s=1.5106 N/m，而混合区域孔隙率为 0.33。探头贯入方法和围压与 McDowell 等 8相同，以 20 mm/s 的速度向下贯入，施加 300 kPa 的围压。模拟结果如图 14，15 所示。将曲线按照下式定义各组模拟结果与 A 组结果的平均残差：， (4)MA2e1()Niiix式中， 分别为某模拟结果曲线上某点的值和Miix和A 组曲线
15、上相对应的值。在曲线均匀取 200 个点进行平均。表 3 CPT 模型几何参数编号 区域数区域外半径/mm区域粒径直径/mm颗粒数量备注A 组 1 150 6 9700B 组 2 75，150 6，12 3550C 组 3 75，95，150 6，9，12 3900文献8D 组 3 75，95，150 6，9，12 3827E 组 3 75，115，150 6，6 + 18，18 3916本文模拟3.1 刚度优化结果对比图 14 为 D 组模拟的端阻应力与深度的结果。可以发现刚度优化后的结果要比大小颗粒采用相同参数的结果更接近与均一小颗粒模拟结果。其中 C 组的平均残差为 0.298，而 D
16、组的平均残差为 0.237。说明，通过将大颗粒区组构与小颗粒区一致，刚度参数基于刚度比和法向刚度与粒径比相同，根据粒径相应优化可有效改善模拟结果，减小颗粒数目。3.2 混合层结果对比图 15 为 E 组模拟的端阻应力与深度的结果。可以发现 B 组由于在变粒径处小颗粒渗入大颗粒中，减弱侧向应力，使得端阻应力要明显小于均一小颗粒模拟结果。通过在变粒径处设置粒径混合层，可有效防止小颗粒渗入，B 组平均残差 0.436，E 组平均残差 0.286，E 组模拟结果显然较 B 组更接近均一小颗粒结果，改善了模拟结果。图 2 颗粒从初始边界的迁移图 14 D 组端阻应力 qc 与深度模拟结果图 15 E 组
17、端阻应力 qc 与深度模拟结果4 结 论本文在已有的变粒径建模方法基础上，针对不同粒径区宏观弹性不一致和界面区大小颗粒互相渗入这两个问题进行优化，通过离散元模拟验证了该方法的有效性。主要得出以下结论： （1）基于宏观弹性一致条件，建立静力条件下几何上等比例放大颗粒体系的颗粒材料参数尺度关系，并满足应力、应变和应变能一致性。（2）通 过 离 散 元 模 拟 验 证 了 不 同 粒 径 的 均 一 颗粒 体系宏观弹性一致情况下的颗粒参数尺度关系。（3）针 对 界 面 区 大 小 颗 粒 互 相 渗 入 ， 设 立 混 合层 ， 提出了两种粒径混合生成特定密实度方法， 并通过离散元模拟验证了相同密实
18、度情况下，两种粒径混合区、混合颗粒与均一颗粒搭接体和均一颗粒区宏观变形特性的一致性。粒径混合层的厚度和接触组构对大尺度模拟结果的影响还有待于进一步研究。参考文献：1 周 健, 崔积弘 , 贾敏才, 等. 静力触探试验的离散元数值模拟研究J. 岩土工程学报, 2007, 29(11): 16041610. 2 JIANG M J, YU H S, HARRIS D. Discrete element modelling of deep penetration in granular soilsJ. International Jouranl for Numerical and Analytica
19、l Methods in Geomechanics, 2006, 30: 335361.3 SHI D D, XUE J F, ZHAO Z Y, et al. A DEM investigation on simple shear behavior of dense granular assembliesJ. Journal of Central South University, 2015, 22(12): 48444855.4 LOMMEN S, SCHOTT D, LODEWIJKS G. DEM speedup: stiffness effects on behavior of bu
20、lk materialJ. Particuology, 2014, 12: 107112.5 BUTLANSKA J, ARROYO M. Homogeneity and symmetry in DEM models of cone penetrationC/ American Institute 4 三峡大学学报（自然科学版） 2017 年4of Physics Conference Proceedings. New York, 2009: 425429.6 CUI L, OSULLIVAN C. An analysis of the triaxial apparatus using a m
21、ixed boundary three-dimensional discrete element modelJ. Gotechnique, 2007, 29(11): 10641610.7 邓益兵, 周 健, 刘文白, 等. 螺旋挤土桩下旋成孔过程的颗粒流数值模拟J. 岩土工程学报, 2011, 33(9): 13911398. 8 MCDOWELL G R, FALAGUSH O, YU H S. A particle refinement method for simulating DEM of cone penetration testing in granular materialsJ
22、. Geotechnique Letters, 2012, 2: 141147.9 CHANG C S, CHAO S J, CHANG Y. Estimates of elastic moduli for granular material with anisotropic random packing structureJ. International Journal of Solids and Structures, 1995, 32(14): 19892008.10 HICHER P Y, CHANG C S. Evaluation of two homogenization tech
23、niques for modeling the elastic behavior of granular materialsJ. Journal of Engineering Mechanics, 2005, 131(11): 11841194.11 ASTM D 4253-00 2002. Standard test method for maximum index density and unit weight of soils using a vibratory tableS. 2002.12 AN X Z, LI C X, YANG R Y, et al. Experimental s
24、tudy of the packing of mono-sized spheres subjected to one-dimensional vibrationJ. Powder Technology, 2009, 196: 5055.13 SHEN S, YU H. Characterized packing of aggregate particles for paving materials: particle size impactJ. Constr Build Mater, 2011, 25(3): 13621368.14 CHO G C, DODDS J, SANTAMARINA
25、J C. Particle shape effects on packing density, stiffness, and strength: natural and crushed sandsJ. Journal of Geotechnical and Geoenvironmental Engineering, 2006, 132(5): 591602.15 MCGEARY R K. Mechanical packing of spherical particlesJ. Journal of the American ceramic Society, 1961, 44(10): 51352
26、2.16 LADE P V, LIGGIO C D, YAMAMURO J A. Effects of non-plastic fines on minimum and maximum void ratios of sandJ. Geotechnical Testing Journal, 1998, 21(4): 336347.17 CHANG C S, MEIDANI M. Dominant grains network and behavior of sandsilt mixtures: stressstrain modelingJ. International Journal for N
27、umerical and Analytical Methods in Geomechanics, 2013, 37: 25632589.18 CHANG C S, WANG J Y, GE L. Modeling of minimum void ratio for sandsilt mixturesJ. Engineering Geology, 2015, 196: 293304.19 CUBRINOVSKI M, ISHIHARA K. Maximum and minimum void ratio characteristics of sandsJ. Soils and Foundations, 2002, 42(6): 6578.
