《三峡大学学报(自然科学版)》：三维离散元大尺度模拟中变粒径方法.doc

第39卷第1期三峡大学学报(自然科学版) 2017年月 Journal of China Three Gorges University(Natural Sciences) DOle 三维离散元大尺度模拟中变粒径方法作者11,作者221,作者31,作者41 (1.上海海事大学海洋科学与工程学院,上海201306:2.中交上海航道勘察设计研究院有限公司,上海200120) 摘要:在已有的变粒径建模方法基础上,针对不同粒径区宏观弹性不一致和界面区大小颗粒互相渗入两个问题进行优化。通过在近场核心区域采用小粒径颗粒模拟真实土颗粒,远场边界区域逐级放大颗粒粒径,不同粒径之间采用相同密实度的混合粒径模拟,实现三维大尺度离散元模拟,以减少模拟颗粒数量提升运算速度。最后,通过模拟静力触探试验(CPT),表明该方法的有效性。该方法可为三维离散元模拟工程大尺度问题提供新途径关键词:离散元模拟:变粒径:弹性常数:尺度关系:相对密实度:静力触探中图分类号:TU43 文献标识码:A 文章编号:1000-4548201701-0062-09 作者简介:作者1(1980-),男,讲师,博士,主要从事土体细观力学与离散元模拟等方面的教学和科研。E-mail: Refinement and application of variable particle-size methods in 3D discrete element modelling for large-scale problems Author l Author 2, Author 3 Author 4 1. College of Ocean and Science Engineering, Shanghai Maritime University, Shanghai 201306, China, 2. Shanghai Waterway Engineering Design and Consulting Co, Ltd. Shanghai 200120, China) 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 sed in the near field core area to simulate 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. For scaling particle systems, the scale relationship is estab lished on the static condition by the identical macro elasticity of uniform granules. This study may provide a ch to the application of 3d discrete element modelling in the large-scale geotechnical engineering. Key discrete element modeling; variable particle-size, elastic constant, scaling law; relative density; CPT 0引言问题划分若干区域,如图1所示,从近场到远场划分砂土具有夭然散粒特性,虽然传统的连续介质力为Z,Z,…,Z区域,在近场核心区域采用小颗粒学方法利用本构方程可以模拟一些工程问题,但对于模拟,以模拟其复杂的物理力学特性;其它各区域模局部的大变形破坏区域和土与结构非连续接触作用区拟颗粒粒径逐级放大,以模拟边界。同时在相邻变粒域,连续性假设不再有效。离散元利用颗粒间的接触径区域增设混合区域M12,M3,…,M(rm,防止相作用可从颗粒尺度上模拟散粒材料的力学行为,目前,邻粒径交界处出现相互嵌入混合。混合区为相邻区域成为研究岩土材料大变形破坏等问题的重要方法。早的两种颗粒粒径混合,密实度与相邻区域一致期学者们多采用二维圆盘离散元模拟13,但二维圆盘由于各区域粒径进行了放大,需要按照尺度关系运动约束不同于真实的颗粒。而三维离散元从颗粒尺调整颗粒材料参数,以保证不同区域相应的物理力学度上模拟工程尺度问题往往需要大量颗粒数目,如何特性的一致性。提高计算效率成为学者们所关注的问题。基金项目:国家自然科学基金青年基金项目() 收稿日期:·通讯作者( E-mail:) 1变粒径建模方法在DEM模拟中,为了减小颗粒数目,将所模拟

第 39 卷第 1 期三峡大学学报（自然科学版） Vol.39 No.1 2017 年 .1 月 Journal of China Three Gorges University (Natural Sciences) Jan. 2017 三维离散元大尺度模拟中变粒径方法作者 11，作者 22，1，作者 3*1，作者 41 (1. 上海海事大学海洋科学与工程学院，上海 201306；2. 中交上海航道勘察设计研究院有限公司，上海 200120) 摘要：在已有的变粒径建模方法基础上，针对不同粒径区宏观弹性不一致和界面区大小颗粒互相渗入两个问题进行优化。通过在近场核心区域采用小粒径颗粒模拟真实土颗粒，远场边界区域逐级放大颗粒粒径，不同粒径之间采用相同密实度的混合粒径模拟，实现三维大尺度离散元模拟，以减少模拟颗粒数量提升运算速度。最后，通过模拟静力触探试验（CPT），表明该方法的有效性。该方法可为三维离散元模拟工程大尺度问题提供新途径。关键词：离散元模拟；变粒径；弹性常数；尺度关系；相对密实度；静力触探中图分类号：TU43 文献标识码：A 文章编号：1000–4548(2017)01–0062–09 作者简介：作者 1(1980– )，男，讲师，博士，主要从事土体细观力学与离散元模拟等方面的教学和科研。E-mail: Refinement and application of variable particle-size methods in 3D discrete element modelling for large-scale problems Author11 , 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, China) 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 simulate 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. 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: discrete element modeling; variable particle-size; elastic constant; scaling law; relative density; CPT 0 引言砂土具有天然散粒特性，虽然传统的连续介质力学方法利用本构方程可以模拟一些工程问题，但对于局部的大变形破坏区域和土与结构非连续接触作用区域，连续性假设不再有效。离散元利用颗粒间的接触作用可从颗粒尺度上模拟散粒材料的力学行为，目前，成为研究岩土材料大变形破坏等问题的重要方法。早期学者们多采用二维圆盘离散元模拟[1-3]，但二维圆盘运动约束不同于真实的颗粒。而三维离散元从颗粒尺度上模拟工程尺度问题往往需要大量颗粒数目，如何提高计算效率成为学者们所关注的问题。 1 变粒径建模方法在 DEM 模拟中，为了减小颗粒数目，将所模拟问题划分若干区域，如图 1 所示，从近场到远场划分为 Z1，Z2，...，Zn 区域，在近场核心区域采用小颗粒模拟，以模拟其复杂的物理力学特性；其它各区域模拟颗粒粒径逐级放大，以模拟边界。同时在相邻变粒径区域增设混合区域 M12，M23，...，M(n-1)n，防止相邻粒径交界处出现相互嵌入混合。混合区为相邻区域的两种颗粒粒径混合，密实度与相邻区域一致。由于各区域粒径进行了放大，需要按照尺度关系调整颗粒材料参数，以保证不同区域相应的物理力学特性的一致性。 ─────── 基金项目：国家自然科学基金青年基金项目（）；收稿日期： *通讯作者（E-mail:） DOI：

峡大学学报(自然科学版) 017年 2不同粒径颗粒体系的弹性验证 8×10 8×105 2.1粒径不同、细观结构相同试样生成方法 8×107 8×10 0.8 4302.4×10824×1050.8同“2 将某一区域小颗粒体系放大至另一区域中,颗粒“5 4.8×1 粒径放大λ倍,按下式构建大颗粒体系坐标系 63024×10824×10508同1 r=x+c 式中,c为坐标平移矢量。 "1,"2,3为随机生成试样:"6,7与”1试样接触组在大颗粒区域生成与小颗粒相对应的大颗粒,如致,法向刚度与粒径比保持相同。图3所示,粒径为小颗粒的λ倍,位置矢量为相应小颗粒位置矢量的λ倍,即 3静力触探变粒径DEM模拟采用上述变粒径方法模拟静力触探,与全部小颗 (2)粒模拟结果对比分析,以验证变粒径方法的有效性式中,x"和x"分别为大和小颗粒体系中第n个球的根据魏杰20理论推导研究,塑性区域在4D(D 位置矢量为探杆直径)范围内。 Huang等的离散元静力触探模拟研究表明探头对土体的应力应变响应主要也集中在4D范围以内。因此在本文的模型中小颗粒集合体区域半径为75mm>4D(72mm),模拟静力触探大变形,75~150mm区域呆用变粒径模拟。模拟分区情况和各区域采用的颗粒粒径列于表3,A组、B组和C组为 McDowel等结果,A组采用均一小颗粒模拟,将其结果作为参照,B组和C组分区分别采用 2种和3种粒径模拟,但大小颗粒细观参数相同,法图1颗粒粒径放大示意图向接触刚度k=5×105Nm,切向接触刚度k=5×105 对于任意相接触的a,b两球,小、大颗粒中心Nm,摩擦系数为0.5,试样孔隙率为03545。D组和连线矢量分别为 E组采用3.1节所述方法生成组构与小颗粒区相同的大颗粒体系,同时,通过保持刚度比和法向刚度与粒 I=x-x=i(x-x=a (3)径比相同,根据颗粒粒径,优化大颗粒刚度参数。D 说明此时大小颗粒体系的接触分布一致,即接触组采用3种粒径模拟细观参数分别为k=5×105N/m, 组构相同 k=5×105N/m:kn=75×105N/m,k=7.5×105N/m和 2.2三轴试验模拟验证 k1=1×106N/m,k=1×10°N/m,E组按4.1节所述方法通过DEM模拟三轴试验,比较分析不同粒径颗在大小颗粒区中间设置一层粒径混合区(见图13), 粒体系的宏观弹性特性,验证不同粒径颗粒体系宏观其颗粒细观刚度为k=5×105Nm,k=5×10Nmn和弹性一致情况下的颗粒参数尺度关系。模拟的工况列 kn=1.5×10N/m,k=1.5×10N/m,而混合区域孔隙于表1中。随机生成3种不同粒径(10,30,60mm)、率为0.33。探头贯入方法和围压与 McDowell等相孔隙比和尺寸相同的三轴试样(见图4)。#1,“2和3同,以20mms的速度向下贯入,施加300kPa的围试验中,3种试样采用相同的细观参数。#4和#5试验压。模拟结果如图14,15所示。将曲线按照下式定义中,保持刚度比α相同,按照法向刚度与粒径比与艹1各组模拟结果与A组结果的平均残差: 中相同,相应增大刚度。“6和试样按照上述粒径不 1、(x-x2)2 同、接触组构相同试样生成方法,在*1试样基础上等例放大,生成与1组构相同试样,按照法向刚度与式中,x和分别为某模拟结果曲线上某点的值和A 粒径比与1中相同,相应增大刚度。施加围压5MPa组曲线上相对应的值。在曲线均匀取200个点进行平进行小应变弹性模拟,采用线弹性接触模拟接触作用模拟结果见图5。表3CPT模型几何参数表1模拟工况编号区域区域外半径区域粒径直径颗粒备注编粒径法向刚度号/mmk/(Nm 切向刚度摩擦组构 k/Nm)系数 A组1

2 三峡大学学报（自然科学版） 2017 年 2 2 不同粒径颗粒体系的弹性验证 2.1 粒径不同、细观结构相同试样生成方法将某一区域小颗粒体系放大至另一区域中，颗粒粒径放大 λ 倍，按下式构建大颗粒体系坐标系： x x c  = + ， (1) 式中， c 为坐标平移矢量。在大颗粒区域生成与小颗粒相对应的大颗粒，如图 3 所示，粒径为小颗粒的  倍，位置矢量为相应小颗粒位置矢量的  倍，即 n n n n d' d  =    =  ， x x ， (2) 式中， n x  和 n x 分别为大和小颗粒体系中第 n 个球的位置矢量。图 1 颗粒粒径放大示意图对于任意相接触的 a，b 两球，小、大颗粒中心连线矢量分别为 b a b a b a ( ) c c c   = −      = − = − =  ，。 l x x l x x x x l (3) 说明此时大小颗粒体系的接触分布一致，即接触组构相同。 2.2 三轴试验模拟验证通过 DEM 模拟三轴试验，比较分析不同粒径颗粒体系的宏观弹性特性，验证不同粒径颗粒体系宏观弹性一致情况下的颗粒参数尺度关系。模拟的工况列于表 1 中。随机生成 3 种不同粒径（10，30，60 mm）、孔隙比和尺寸相同的三轴试样（见图 4）。#1，#2 和#3 试验中，3 种试样采用相同的细观参数。#4 和#5 试验中，保持刚度比  相同，按照法向刚度与粒径比与#1 中相同，相应增大刚度。#6 和#7 试样按照上述粒径不同、接触组构相同试样生成方法，在#1 试样基础上等比例放大，生成与#1 组构相同试样，按照法向刚度与粒径比与#1 中相同，相应增大刚度。施加围压 5 MPa 进行小应变弹性模拟，采用线弹性接触模拟接触作用，模拟结果见图 5。表 1 模拟工况编号粒径 /mm 法向刚度 kn/(N·m -1 ) 切向刚度 ks//(N·m -1 ) 摩擦系数组构 # 1 10 8×107 8×105 0.8 — # 2 30 8×107 8×105 0.8 — # 3 60 8×107 8×105 0.8 — # 4 30 2.4×108 2.4×106 0.8 同# 2 # 5 60 4.8×108 4.8×106 0.8 同# 3 # 6 30 2.4×108 2.4×10 6 0.8 同# 1 # 7 60 4.8×108 4.8×106 0.8 同# 1 注：# 1，# 2，# 3 为随机生成试样；# 6，# 7 与# 1 试样接触组构一致，法向刚度与粒径比保持相同。 3 静力触探变粒径 DEM 模拟采用上述变粒径方法模拟静力触探，与全部小颗粒模拟结果对比分析，以验证变粒径方法的有效性。根据魏杰[20]理论推导研究，塑性区域在 4D（D 为探杆直径）范围内。Huang 等[21]的离散元静力触探模拟研究表明探头对土体的应力应变响应主要也集中在 4D 范围以内。因此在本文的模型中小颗粒集合体区域半径为 75 mm>4D（72 mm），模拟静力触探大变形，75～150 mm 区域采用变粒径模拟。模拟分区情况和各区域采用的颗粒粒径列于表 3，A 组、B 组和 C 组为 McDowell 等[8]结果，A 组采用均一小颗粒模拟，将其结果作为参照，B 组和 C 组分区分别采用 2 种和 3 种粒径模拟，但大小颗粒细观参数相同，法向接触刚度 kn=5×105 N/m，切向接触刚度 ks=5×105 N/m，摩擦系数为 0.5，试样孔隙率为 0.3545。D 组和 E 组采用 3.1 节所述方法生成组构与小颗粒区相同的大颗粒体系，同时，通过保持刚度比和法向刚度与粒径比相同，根据颗粒粒径，优化大颗粒刚度参数。D 组采用 3 种粒径模拟细观参数分别为 kn=5×105 N/m， ks=5×105 N/m；kn=7.5×105 N/m，ks=7.5×105 N/m 和 kn=1106 N/m，ks=1×106 N/m，E 组按 4.1 节所述方法在大小颗粒区中间设置一层粒径混合区（见图 13），其颗粒细观刚度为 kn=5×105 N/m，ks=5×105 N/m 和 kn=1.5×106 N/m，ks=1.5×106 N/m，而混合区域孔隙率为 0.33。探头贯入方法和围压与 McDowell 等[8]相同，以 20 mm/s 的速度向下贯入，施加 300 kPa 的围压。模拟结果如图 14，15 所示。将曲线按照下式定义各组模拟结果与 A 组结果的平均残差： M A 2 e A 2 1 1 ( ) ( ) N i i i i x x N x  = − =  ， (4) 式中， M A i i x x 和分别为某模拟结果曲线上某点的值和 A 组曲线上相对应的值。在曲线均匀取 200 个点进行平均。表 3 CPT 模型几何参数编号区域数区域外半径 /mm 区域粒径直径 /mm 颗粒数量备注 A 组 1 150 6 9700

组275,1506,123550文献8]4结论 C组375,95,1506,9,123900 本文在已有的变粒径建模方法基础上,针对不同 D组375,95,1506,9,123827本文模粒径区宏观弹性不一致和界面区大小颗粒互相渗入这 E组375,115,1506,6+18,183916抄两个问题进行优化,通过离散元模拟验证了该方法的 3.1刚度优化结果对比有效性。主要得出以下结论图14为D组模拟的端阻应力与深度的结果。可 (1)基于宏观弹性一致条件,建立静力条件下几以发现刚度优化后的结果要比大小颗粒采用相同参数何上等比例放大颗粒体系的颗粒材料参数尺度关系的结果更接近与均一小颗粒模拟结果。其中C组的平并满足应力、应变和应变能一致性。均残差为0.298,而D组的平均残差为0237。说明, (2)通过离散元模拟验证了不同粒径的均一颗粒通过将大颗粒区组构与小颗粒区一致,刚度参数基于体系宏观弹性一致情况下的颗粒参数尺度关系刚度比和法向刚度与粒径比相同,根据粒径相应优化 (3)针对界面区大小颗粒互相渗入,设立混合层, 可有效改善模拟结果,减小颗粒数目。提出了两种粒径混合生成特定密实度方法,并通过离 32混合层结果对比散元模拟验证了相同密实度情况下,两种粒径混合区、图15为E组模拟的端阻应力与深度的结果。可混合颗粒与均一颗粒搭接体和均一颗粒区宏观变形特以发现B组由于在变粒径处小颗粒渗入大颗粒中,减性的一致性弱侧向应力,使得端阻应力要明显小于均一小颗粒模粒径混合层的厚度和接触组构对大尺度模拟结果拟结果。通过在变粒径处设置粒径混合层,可有效防的影响还有待于进一步研究止小颗粒渗入,B组平均残差0.436,E组平均残差 0.286,E组模拟结果显然较B组更接近均一小颗粒结参考文献: 果,改善了模拟结果。 [周健,崔积弘,贾敏才,等.静力触探试验的离散元数值模拟研究[岩土工程学报,2007,29(11):1604-1610. 2JJIANG MJ, YUHS, HARRIS D Discrete element modelling or Numerical and Analytical Methods in Geomechanics, 图2颗粒从初始边界的迁移 2006,30:335-361 0 [3]SHID D, XUE J F, ZHAO ZY, et al. A DEM investigation on simple shear behavior of dense granular assemblies] ty,2015,22(12) 4 LOMMEN S, SCHOTT D, LODEWIJKS G. DEM speedup stiffness effects on behavior of bulk material[]. Particuology, 2014.12:107-112. 5]BUTLANSKAJ, ARROYOM. Homogeneity and symmetry in DEM models of cone penetration[ C]/ American Institute of 图14D组端阻应力φ与深度模拟结果 Physics Conference Proceed ings. New York. 2009: 425 6 CUI L, O'SULLIVAN C. An analysis of the triaxial apparatus using a mixed bound ary three-dimensional discrete element modelp]. Geotechnique, 2007, 29(11): 1064-1610 7]邓益兵,周健,刘文白,等.螺旋挤土桩下旋成孔过程的颗粒流数值模拟[.岩土工程学报,2011,33(9):1391 8 MCDOWELL G R, FALAGUSH O, YU H S. A particle 9o/MPa efinement method for simulating DEM of cone penetration 图15E组端阻应力φ与深度模拟结果 testing in granular mater ials]. Geotechnique Letters, 2012

B 组 2 75，150 6，12 3550 文献[8] C 组 3 75，95，150 6，9，12 3900 D 组 3 75，95，150 6，9，12 3827 本文模 E 组 3 75，115，150 6，6 + 18，18 3916 拟 3.1 刚度优化结果对比图 14 为 D 组模拟的端阻应力与深度的结果。可以发现刚度优化后的结果要比大小颗粒采用相同参数的结果更接近与均一小颗粒模拟结果。其中 C 组的平均残差为 0.298，而 D 组的平均残差为 0.237。说明，通过将大颗粒区组构与小颗粒区一致，刚度参数基于刚度比和法向刚度与粒径比相同，根据粒径相应优化可有效改善模拟结果，减小颗粒数目。 3.2 混合层结果对比图 15 为 E 组模拟的端阻应力与深度的结果。可以发现 B 组由于在变粒径处小颗粒渗入大颗粒中，减弱侧向应力，使得端阻应力要明显小于均一小颗粒模拟结果。通过在变粒径处设置粒径混合层，可有效防止小颗粒渗入，B 组平均残差 0.436，E 组平均残差 0.286，E 组模拟结果显然较 B 组更接近均一小颗粒结果，改善了模拟结果。图 2 颗粒从初始边界的迁移图 14 D 组端阻应力 qc 与深度模拟结果图 15 E 组端阻应力 qc 与深度模拟结果 4 结论本文在已有的变粒径建模方法基础上，针对不同粒径区宏观弹性不一致和界面区大小颗粒互相渗入这两个问题进行优化，通过离散元模拟验证了该方法的有效性。主要得出以下结论：（1）基于宏观弹性一致条件，建立静力条件下几何上等比例放大颗粒体系的颗粒材料参数尺度关系，并满足应力、应变和应变能一致性。（2）通过离散元模拟验证了不同粒径的均一颗粒体系宏观弹性一致情况下的颗粒参数尺度关系。（3）针对界面区大小颗粒互相渗入，设立混合层，提出了两种粒径混合生成特定密实度方法，并通过离散元模拟验证了相同密实度情况下，两种粒径混合区、混合颗粒与均一颗粒搭接体和均一颗粒区宏观变形特性的一致性。粒径混合层的厚度和接触组构对大尺度模拟结果的影响还有待于进一步研究。参考文献： [1] 周健, 崔积弘, 贾敏才, 等. 静力触探试验的离散元数值模拟研究[J]. 岩土工程学报, 2007, 29(11): 1604–1610. [2] JIANG M J, YU H S, HARRIS D. Discrete element modelling of deep penetration in granular soils[J]. International Jouranl for Numerical and Analytical Methods in Geomechanics, 2006, 30: 335–361. [3] SHI D D, XUE J F, ZHAO Z Y, et al. A DEM investigation on simple shear behavior of dense granular assemblies[J]. Journal of Central South University, 2015, 22(12): 4844– 4855. [4] LOMMEN S, SCHOTT D, LODEWIJKS G. DEM speedup: stiffness effects on behavior of bulk material[J]. Particuology, 2014, 12: 107–112. [5] BUTLANSKA J, ARROYO M. Homogeneity and symmetry in DEM models of cone penetration[C]// American Institute of Physics Conference Proceedings. New York, 2009: 425– 429. [6] CUI L, O’SULLIVAN C. An analysis of the triaxial apparatus using a mixed boundary three-dimensional discrete element model[J]. Géotechnique, 2007, 29(11): 1064–1610. [7] 邓益兵, 周健, 刘文白, 等. 螺旋挤土桩下旋成孔过程的颗粒流数值模拟[J]. 岩土工程学报, 2011, 33(9): 1391– 1398. [8] MCDOWELL G R, FALAGUSH O, YU H S. A particle refinement method for simulating DEM of cone penetration testing in granular materials[J]. Geotechnique Letters, 2012