（岩土工程专业论文）基于coulomb摩擦阻尼的离散单元法.pdf

西南交通大学硕士研究生学位论文第1 页 摘要 土是由大小不同、形状各异的颗粒构成的不连续介质。作为非连续介质的计算方 法，离散单元法非常适合求解颗粒介质问题。该法中，通常采用粘性阻尼耗散系统运 动的能量，从而使计算快速收敛。但实际上，颗粒之间的切向作用力更接近c o u l o m b 摩擦力。因此，采用基于c o u l o m b 摩擦阻尼( 即摩擦系数) 建立的离散单元法接触模 型，能够真实地描述颗粒间的切向作用，并得到整个系统真实的稳定时间。本文将尝 试对现有离散单元法计算模型进行改进，建立基于c o u l o m b 摩擦阻尼的接触模型及相 应的计算方法。 论文首先计算分析了一维振动问题，给出滑块在静止或滑动状态下摩擦力的确定 方法及运动状态的判定方法。并通过1 个滑块和5 0 个滑块在接触面上滑动的例子，说 明了基于c o u l o m b 摩擦阻尼运动过程的数值法求解原理及方法。对二维问题，建立基 于c o u l o m b 摩擦阻尼的接触模型，推导了相应的离散单元法的计算公式。其中，重点 解决了颗粒间切向摩擦力的计算方法。上述工作完成后，编制了相应的计算程序，并 计算了一个简单的算例，验证了本文计算模型的合理性及计算方法的可行性。 与传统的离散单元法相比，基于c o u l o m b 摩擦阻尼的离散单元法具有如下优点： ( 1 ) c o u l o m b 摩擦更为真实地描述了颗粒之间在接触点的切向作用。( 2 ) 以c o u l o m b 摩擦阻尼代替粘性阻尼，其值更易确定，且由此得到的稳定时间是真实的时间。( 3 ) 只有颗粒间的法向刚度系数而没有切向刚度系数，计算参数的确定更为简单。 关键词：离散单元法；颗粒；c o u l o m b 摩擦阻尼；切向摩擦力 西南交通大学硕士研究生学位论文第1 | 页 曼曼皇曼曼曼曼曼皇曼! 曼曼曼曼曼曼皇曼鼍皇曼鼍皇曼曼i i m m 一! m m l 鼍曼曼皇曼曼鼍皇曼曼曼! 曼曼! 曼曼曼曼皇曼曼曼曼曼鼍曼曼皇曼曼皇曼曼 a b s t r a c t s o i li sd i s c o n t i n u o u sm e d i u mw h i c hi sc o m p o s e do fg r a n u l e sw i t hd i f f e r e n ts i z e sa n d v a r i e ds h a p e a san u m e r i c a lm e t h o df o rd i s c o n t i n u o u sm e d i u m ，d i s c r e t ee l e m e n tm e t h o di s p e r f e c t l ys u i t a b l et os o l v et h ep r o b l e m sr e l a t e dt o t h eg r a n u l a rm e d i a u s u a l l y , v i s c o u s d a m p i n gi su s e db yd e m t od i s s i p a t et h ee n e r g yi nt h em o v e m e n to fg r a n u l a rs y s t e ms oa s t of o r c et h ec a l c u l a t i o nt ob ec o n v e r g e dq u i c k l y a sam a t t e ro ff a c t ，t a n g e n t i a lc o n t a c tf o r c e b e t w e e ng r a n u l e sa c c o r d sw i t hc o u l o m bf r i c t i o nf o r c e i nt h i sp a p e r , ac o n t a c tm o d e lb a s e d o nc o u l o m bd a m p i n gw a sa p p l i e di nd i s c r e t ee l e m e n tm e t h o di n s t e a do fv i s c o u sd a m p i n g ， a n dr e l e v a n tc o m p u t a t i o n a lm e t h o d sw i t l lc o u l o m bd a m p i n gw a sd e v e l o p e d f i r s t l y , t h ec h a r a c t e r i s t i c so fo n e - d i m e n s i o n a lv i b r a t i o nw i t hc o u l o m bd a m p i n gw e r e a n a l y z e d t h en u m e r i c a lm e t h o dt os o l v ev i b r a t i o nw i t hc o u l o m bd a m p i n gw a si l l u s t r a t e d w i mt h ev i b r a t i o no fab l o c ka n dag r o u po fb l o c k so naf r i c t i o n a ls u r f a c e a sf o rt h e t w o - d i m e n t i o a lp r o b l e m s ，t h ec o n t a c tm o d e lo nb a s i so fc o u l o m bd a m p i n gw a se s t a b l i s h e d a n dt h ec o r r e s p o n d i n gc o m p u t a t i o n a lf o r m u l ao fd i s c r e t ee l e m e n tm e t h o dw a sd e r i v e d t h e c a l c u l a t i o n a lm e t h o d so ft h et a n g e n t i a lf r i c t i o nb e t w e e ng r a n u l e sw e r eo n eo fd i f f i c u l t p r o b l e m si nt h i sm e t h o d t h ep r o g r a mw a sc o m p i l e dw i t hc * l a n g u a g e as i m p l ee x a m p l e w a sc o m p u t e dt oi l l u s t r a t et h et h er a t i o n a l i t yo fc a l c u l a t i o nm o d e l sa n df e a s i b i l i t yo f c a l c u l a t i o n a lm e t h o d si nt h i sp a p e r c o m p a r e dw i t ht h et r a d i t i o n a ld e m ，d i s c r e t ee l e m e n tm e t h o db a s e do nc o u l o m b d a m p i n gh a st h ef o l l o w i n ga d v a n t a g e s f i r s t l y , c o u l o m bd a m p i n gi sa b l et om o r et r u e l y d e s c r i b et h et a n g e n t i a le f f e c t sb e t w e e ng r a n u l e s s e c o n d l y , t h em a t e r i a lc o e f f i c i e n t si nt h e c o u l o m bd a m p i n gm o d e li se a s i e rt ob ed e t e r m i n e dt h a nt h o s ei nc o n v e n t i o n a ld e m m o r e o v e r , t h e r ei so n l yc o e f f i c i e n to fn o r m a ls t i f f n e s sr a t h e rt h a nc o e f f i c i e n to ft a n g e n t i a l s t i f f n e s sb e t w e e ng r a n u l e si nt h em o d e l ，w h i c hi sr a t h e rc o n v e n i e n tf o rd e t e r m i n a t i o no f c o m p u t a t i o n a lp a r a m e t e r s k e yw o r d s ：d i s c r e t ee l e m e n tm e t h o d ；g r a n u l e s ；c o u l o m bd a m p i n g ；t a n g e n t i a l f r i c t i o nf o r c e 西南交通大学曲南父逋大罕 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保留并 向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅。本人授 权西南交通大学可以将本论文的全部或部分内容编入有关数据库进行检索，可以采用 影印、缩印或扫描等复印手段保存和汇编本学位论文。 本学位论文属于 1 保密口，苟 年解密后适用本授权书； 2 不保密彤使用本授权书。 ( 请在以上方框内打“ ) 学位论文作者签名：匆辟辱臣 日期：知f d i ( ，。) 西南交通大学硕士学位论文主要工作( 贡献) 声明 本人在学位论文中所做的主要工作或贡献如下： ( 1 ) 对现有离散单元法的计算模型进行改进，以c o u l o m b 摩擦阻尼替代粘性阻尼， 用以描述颗粒间切向作用，更为真实地反映了颗粒介质的受力特点。同时，计算模型 中不需引入切向刚度，计算参数的确定更为简单。 ( 2 ) 基于c o u l o m b 摩擦阻尼模型，建立了二维离散单元法的计算公式及方法，并 编制了相应的计算程序。 本人郑重声明：所呈交的学位论文，是在导师指导下独立进行研究工作所得的成 果。除文中已经注明引用的内容外，本论文不包含任何其他个人或集体已经发表或撰 写过的研究成果。对本文的研究做出贡献的个人和集体，均已在文中作了明确说明。 本人完全了解违反上述声明所引起的一切法律责任将由本人承担。 学位论文作者签名：硼霉臣 日期：沏m 6 7 西南交通大学硕士研究生学位论文第1 页 曼皇曼曼曼皇! 曼曼曼! 曼! 曼曼! 曼i i ml_ i ：a m ： _ _ i 一一i i 篡 第1 章绪论 1 1 研究背景及意义 土是由颗粒组成的不连续体，但在目前的计算和分析中，多将土体作为连续体来 计算分析。针对连续介质的数值分析方法有限差分法、有限元法和边界元法等，上述 计算方法目前已经发展得比较成熟，并被广泛地应用在岩土工程领域。而创建于上世 纪7 0 年代的离散单元法则是以颗粒、块体介质为研究对象的一种数值计算方法，经过 了3 0 多年的发展，该方法已经发展得相当成熟。在理论上，从单纯的离散元模拟发展 到与有限元和边界元的耦合，从二维计算发展到三维计算，从静力学问题的模拟发展 到对动态问题的模拟，从纯固体介质模拟发展到对多相介质、多场耦合问题的模拟等。 在离散单元法中，颗粒( 块体) 之间接触作用是计算模型的关键问题。现有的计 算方法中，采用弹簧阻尼器模型，即：用法( 切) 向接触力+ 粘性阻尼来模拟颗粒( 块体) 之间的作用，法向应力、切向应力分别由法( 切) 向刚度系数和变形确定。而粘性阻尼是 与速度呈正比关系，它主要被用来耗散系统的动能，从而能使计算快速收敛。显然， 这里的阻力只是一个计算参数，它所反映的并不是颗粒之间发生相互位移时所受到的 真正的阻力特性，由此计算模型得到的颗粒介质的最终稳定时间也是一个虚拟量。 事实上，颗粒之间在受力变形过程中，切向的作用力更接近于c o u l o m b 摩擦阻尼， 其特点是阻尼与颗粒之间法向接触力呈正比关系，方向取决于颗粒间的相对运动，这 与粘性阻尼明显不同。而且，从理论上讲，基于c o u l o m b 摩擦阻尼得到的最终稳定时 间是颗粒介质真正的稳定时间。 本论文针对颗粒介质材料，尝试以c o u l o m b 摩擦阻尼代替现有接触模型中的粘性 阻尼。运用离散单元法的计算原理，建立基于c o u l o m b 摩擦阻尼的离散单元计算方法， 推导新模型下的计算公式。然后利用卅语言编制计算程序，并用算例进行计算，以验 证本文计算模型及方法的合理性和可行性。 1 2 颗粒介质散体力学数值计算方法研究现状 目前，颗粒介质散体力学计算数值方法主要有【1 】：离散元法( d i s c r e t ee l e m e n t m e t h o d ) 、不连续变形分析( d e f o r m e dd i s c r e t ea n a l y s i s ) 方法、光滑粒子法( s m o o t h e d 西南交通大学硕士研究生学位论文第2 页 p a r t i c l eh y d r o d y n a m i c sm e t h o d ) 及改进的光滑粒子法( c s p h ) 、颗粒流法等。 如前所述，离散元法由c u n d a l l 于1 9 7 1 年提出。它的基本思想是把不连续体分离 为刚性单元的集合，并使各个刚性单元满足运动方程，用时步迭代的方法求解各刚性 单元的运动方程，继而求得不连续体的整体运动形态，该方法已在岩土工程等领域得 到广泛应用。 不连续变形分析方法是石根华提出的分析块体系统运动和变形的一种数值方法。 它将岩体按节理面及断层因势划分为不同的块体单元，以各单元体的位移为未知量， 根据块体的边界条件及颗粒接触关系建立总体平衡方程，该方程可以由系统的最小势 能原理求得。通过求解方程组即可得到块体在当前时步下的位移、受力及变形等物理 量。通过循环计算求解方程，就可以得到系统达到静力平衡时的位移场及应力场等。 此方法被广泛应用于求解岩土工程中的边坡及滑坡稳定性问题，效果明显。 光滑粒子法( s m o o t h e dp a r t i c l eh y d r o d y n a m i c s ) 是一种典型的l a g r a n g e 无网格算法。 它通过带质量的粒子离散计算域，粒子即代表材料，不同的材料的粒子即构成界面， 不同材料的粒子相对滑动就形成所谓的界面滑移。通过引入了核函数( k e r n e lf u n c t i o n ) 来构造局部光滑的连续场，s p h 将离散体系转化为准连续体系。核函数的意义是用来 表征节点及其影响域( i n f l u e n c ed o m a i n s u p p o r td o m a i n ) l 为物理量间的关系，提供了一种 插值思路。通过直接对强形式的偏微分方程求解，得到系统的运动形态。s p h 在求解 爆炸冲击及大变形问题方面有不少应用【2 1 。 颗粒流法指的是颗粒材料在外力作用和内部应力状况变化时发生的类似于流体的 运动状态。其主要理论是把敖粒体看作有限离散颗粒单元的集合体，将材料理想化为 相互独立、相互接触和相互作用的颗粒群体。对于固相颗粒通过求解运动和动量方程 模拟颗粒运动，对于液相介质采用均一化流体c f d 计算技术模拟其在孔隙中的运动， 即通过求解平均n a v i e r - s t o k e s 方程模拟孔隙中流体的运动。该理论研究始于2 0 世纪 5 0 年代，早先的工作主要集中于一些简单的试验观测和近似分析，并为工程设计提供 一些经验公式。7 0 年代以后，人们开始对颗粒流动的物理机制发生了兴趣，颗粒流研 究进入一个蓬勃发展的阶段【3 1 。 1 3 本文的主要内容及技术路线 如前所述，离散单元法是求解颗粒介质力学问题的一种有效的数值计算方法，但 西南交通大学硕士研究生学位论文第3 页 其描述颗粒接触的力学模型需要改进，所以本论文将以离散单元法的基本理论及计算 方法为基础，尝试将c o u l o m b 摩擦阻尼代替粘性阻尼，建立基于c o u l o m b 摩擦阻尼的 接触模型及计算方法。 相应的技术路线是： ( 1 ) 研究采用c o u l o m b 摩擦阻尼时颗粒相互作用的维问题。分析一维c o u l o m b 摩擦振动的运动特点，为二维问题的实现奠定基础。 ( 2 ) 针对二维问题，建立颗粒间相互作用的力学模型：颗粒简化为圆形，且不发 生转动；颗粒之间法向力与法向相对位移成正比，切向力按c o u l o m b 摩擦力确定。 ( 3 ) 采用与离散单元法相似的方法建立新模型的计算公式。 ( 4 ) 运用卅语言编制相应的计算程序，并选取合适的算例进行计算，对论文计 算模型及计算方法的合理性、正确性进行验证。 西南交通大学硕士研究生学位论文第4 页 曼,h,ii=_ii- i i 鼍曼曼鼍 第2 章离散单元法的计算原理 由于本文将以离散单元法的基本理论为基础，所以本章主要介绍离散单元法的基 本原理。 离散单元法( d i s c r e t ee l e m e n tm e t h o d ) 是c u n d a l l 于2 0 世纪7 0 年代初提出的， 最初它的研究对象主要是岩石等非连续介质的力学行为，将不连续体分离为刚性元素 的集合，使各个刚性元素满足运动方程，用时步迭代的方法求解各刚性元素的运动方 程，继而求得不连续体的整体运动形态，后来在众多学者的努力下，该方法被广泛应 用于求其他非连续介质问题。离散单元法的重要特点是允许单元间的相对运动，甚至 分离，因此非常适合于计算以块体、颗粒等基本单元的散体介质。经过3 0 多年的发展， 离散单元法的计算理论及方法已相当成熟，并在实际工程中得到了广泛的应用。 2 1 离散单元法的发展与应用 c u n d a u 于1 9 7 1 年首次提出了用于分析岩石块体的弹簧阻尼器的接触模型及相应 的计算方法【4 】，又在1 9 7 9 年和s t a c k 5 提出了适合于土力学的软球模型将颗粒间的 接触过程简化为弹簧振子的阻尼振动的模型，并推出了基于二维圆盘的b a l l 程序和 基于三维圆球的t r u b a l 程序【6 1 。1 9 8 0 年，c a m p b e l l 7 1 提出了硬颗粒模型认为在 碰撞时不需要考虑颗粒的变形、接触力等细节，只需要确定碰撞后的速度的模型。1 9 8 9 年，t h o r n t o n 对t r u b a l 进行了全面的改造，抛弃c u n d a l l 的弹簧一阻尼器模型，提 出球体接触力学模型，最终构建了一个新的离散元仿真程序，并编写了新的离散元仿 真程序，命名为g r a n u l e 。这个程序能够模拟干湿、弹性、塑性和颗粒两相流问题， 大大拓宽了离散单元法的应用的范围。1 9 9 9 年，o d a 等 8 】共同编写了专著，较为系统 地介绍了散体力学和离散元法，并对散体宏观细观力学关系的作了研究。 离散元在我国起步较晚，但是发展迅速。1 9 8 6 年，东北大学的王泳嘉 9 首次向我 国岩石力学与工程界介绍了离散元法的基本原理及应用，并在其书中详细介绍了 c u n d a l l 的接触模型，给出了具体的求解过程及关键技术处理，并将该理论应用于巷道 支护、边坡工程、渗流问题等领域。阐述了离散元法与拉格朗日法在岩土力学中的应 用，后来开发了离散元法软件系统2 d b l o c k 、三维离散单元软件t r u d e c 【l 0 1 1 】。 随着离散单元法的不断发展，有学者将其同其他算法融合在一起，其中就包括同 西南交通大学硕士研究生学位论文第5 页 曼! ! ! 曼! ! ! ! ! ! ! ! ! 笪! ! 曼曼! 蔓! ! 曼蔓曼! 曼曼曼曼曼皇皇曼皇曼i ! 曼曼! 曼寰曼皂曼皇曼曼曼曼曼! 曼! 曼! ! 曼曼曼曼曼 有限元、边界元法的耦合，能充分发挥各自的长处，扩大各自算法的应用领域【14 1 。例 如，周晓青等【1 5 】以地下洞室开挖为例讨论了离散元与边界元的耦合计算，认为离散元 和边界元的耦合算法是模拟含节理、裂隙岩体问题的一种有效方法，提出了用离散元 法模拟局部的岩块滑落变形的复杂区域，而以边界元模拟问题域中的无限或半无限弹 性体的方法，提高计算精度，并且减少了计算工作量。 2 1 1 离散单元法在岩土工程中的应用 麻风海等 16 】将上覆岩体视为含有不连续面( 如节理等) 的层状块体模型，采用离散元 法研究了岩层内部直至地表的动态移动过程，并用实测资料得到上覆岩层从直接顶板 直至地表整个范围内的移动动态关系，验证了离散元法的分析结果。谭云亮等【1 7 】针对 巷道围岩多被节理切割的特点，利用离散元法研究了全长锚杆对节理围岩稳定性的影 响，并指出巷道各部位围岩的变形与破坏具有明显协同性或关联性，应予以系统地考 虑对围岩的加固。焦玉勇等【l8 1 9 】利用三维离散元法分析模拟了山地滑坡、地下水和锚 杆，给出了相应的算法和公式。姜卫利等 2 0 】把离散单元法应用于散粒道床的数值分析， 建立道床离散元模型，计算散粒道床在不同加载情况下的应力并对结果进行分析，分 析中还把道渣假设为可变形单元，并在刚性块体模型的基础上，将块体分解成几个常 应变的三角形，采用高斯定理和w i l s o n 的大变形有限差分求解三角形单元的位移和应 变。周立新等【2 l 】用离散元法模拟了岩体卸荷松动模式，并以黄河黑山峡大柳树坝址为 例，概化了坝址区岩体卸荷松动变形模式，同时建立了地质模型。 鲁军【2 2 】在他的三维刚性体离散元法数值模拟中提出了凸多面块体的角边修圆模式 以避免岩块尖锐棱角给计算带来的困难，这对解决不规则形状的颗粒问题将有一定的 启发作用。魏群【2 3 】提出了利用计算机模拟产生任意形状离散颗粒组合体的最新研究成 果，利用蒙特卡诺法，产生满足不同分布函数，不同级别的颗粒组合体，为散体元法 颗粒组合体的变形特性研究提供了基本条件。史旦达【2 4 】等针对纯圆颗粒模拟的缺陷， 利用离散单元法开发了似椭圆的“椭圆团颗粒 ，研究了砂土单调剪切特性及其细观机 理，并与实验结果对比，得到了颗粒形状对试样的强度和变形特性具有重要影响的结 论。吴剑【2 5 】以散粒体力学为基础，研究滑带土的颗粒离散性对其宏观力学性能，尤其 是宏观动力行为的影响，模拟了具有分相特性的滑带离散颗粒集的剪切动力过程，揭 示了具有分相特性滑带在剪切条件下分相差速运动的变化规律。 西南交通大学硕士研究生学位论文第6 页 2 2 离散单元法基本原理 以颗粒介质为例，在离散单元法中，每个颗粒是一个单元，根据颗粒介质组成的 不同，颗粒可以是规则的，也可以是不规则的，理想化的单元为圆球。在计算时，每 一时步扫描颗粒所处位置，根据检索技术，搜索颗粒与周围颗粒的接触情况，当颗粒 间发生接触时，相互之间的作用力可根据由切向和法向的接触模型( 即理论上的弹、塑 性变形量和接触力间的关系) 得出，再根据牛顿运动定律确定颗粒的速度、加速度等物 理量，利用中心差分的方法求解单元在下一时步的位移，如果对所有的单元进行循环 迭代计算，就可以得到系统在下一时步的运动形态。 2 2 1 计算模型及假定 y 图2 1 颗粒的集合及作用于单个颗粒上的力 离散单元法基本假设如下：单元是剐性的，即单元的几何形状不会因单元间的 挤压力作用而改变；由于计算时步间隔取得足够小，单元的速度和加速度在一个时 步内为常量，并且单元在一个时步内只能以很小的位移与其相邻单元作用，其作用力 也只能传递到其邻接单元，而不能传递得更远；单元问的连接是靠相互接触实现的， 圆形单元的接触为点接触。 在离散单元法的模型中，介质被假定为离散块体的集合( 图2 1 ( a ) ) ，假如对某 个颗粒b ，其周围有邻接颗粒通过接触作用于它的一组力( 见图2 - 1 ( b ) ) l ，( f = 1 5 ) 。上述一组力会对颗粒b 产生合力及合力矩，如果合力及合力矩不为 零，则根据牛顿第二定律开始运动，颗粒将发生运动。按照时步迭代方法遍历整个颗 西南交通大学硕士研究生学位论文第7 页 粒集合，直至每个颗粒均不出现不平衡力及1 i 半衡力矩。 2 2 2 离散单元法的基本方程 2 2 2 1 物理方程一颗粒接触力和位移之间关系 如图2 - 2 ( a ) 所示，假定颗粒之间的法向作用力e 正比于其法向“叠合”r u 。，即 只= 吒材。 ( 2 - 1 ) 式中吒法向刚度系数。对于圆形颗粒，法向刚度系数可根据接触颗粒的弹性 模量及接触半径确定； 。接触颗粒之间的法向相对位移。 新位置 ( a ) ( b ) 图2 2 颗粒之间的作用力 如果两个颗粒间发生相对滑移u 。( 见图2 - 2 ( b ) ) ，则接触面上会产生剪切力，即 z x f = k , u ， ( 2 2 ) 式中 颤剪切刚度系数。可根据法向刚度折减得到； u s 颗粒间切向相对滑移值。 2 2 2 2 运动方程 根据牛顿第二定律，对某个颗粒，其加速度为 位移 西南交通大学硕士研究生学位论文第8 页 式中e 、颗粒x 、y 方向的合力； 竹一颗粒的质量。 将式( 2 3 ) 的加速度用向前差分格式表示，可以得到颗粒质心x ，y 方向的速度和 式中a x = i l ； 时步。 2 2 3 求解过程 f v ? 址他一- - k t + 口， 卜芦佗= 1 ，；+ 口，a t i “件& ，= “j + 西，出 弦“y = 甜y + 西y f ( 2 - 4 ) ( 2 5 ) r + a t l 2 时刻速度 t + 删a t 2 t - 州i 2 ：誊篇 任6 ， i “x 2 “x+ 甜j f ，一 、 1 f + 址2 f - 出2 、o o ， 【甜y2 甜y + 甜y f 再对t + a t 2 时刻的速度中心差分后，可得到f + a t 的位移 甜：t+at：甜：t，+：t，+a出t2-jf t(2-7ut + a r 2 a t ) 工工一j ， 、 1r ，+ ， ， l 甜y2 甜y + “y 整理后 ：t+at：“：t，+；西。ttl t + t ytat+西ft：jt一-缸t坨2)出at c 2 - 8 ， i 甜y = “y + ( 西y 、。7 上式中的位移q 及“y 、加速度和呶卜酬2 及舀；和甜 y l - - 驯2 在时步，一a t - t 的计算中 已经确定。 。 擎擎 纵 乃 西南交通大学硕士研究生学位论文第9 页 ( 1 ) 以位移作为基本未知量，根据力与位移关系即式( 2 1 ) 和( 2 2 ) 可以求得 块体间的接触力。 ( 2 ) 由式( 2 3 ) 可以得到块体当前时刻的加速度。 ( 3 ) 由式( 2 6 ) 和( 2 7 ) 进而求得下一时刻的位移、速度等其他物理量，最终 得到整个颗粒系统的运动状态。 具体的求解过程如图2 3 所示，计算路径为一一一一一。在一个时 步中，计算从t 时刻的位移开始，结束于f + 址时刻的位移。 图2 3 离散单元法求解过程示意图 当运动从，= o 时刻开始时，由于一出2 时刻是不存在的，可以使该时刻的速度等 于户o 时刻的速度，然后再根据图2 3 所示的求解方法进行计算。 在解决实际问题时，需要根据边界类型的不同对边界处的颗粒加以限制，常用的 边界控制方法主要有：( 1 ) 刚性壁边界，假设边界是一个刚性墙壁，颗粒位移被限制 在刚性壁以内；( 2 ) 颗粒边界，边界是由颗粒组成的，控制方式有以下4 种：位移控 制边界、应力控制边界、伺服控制边界、边界平均应力控制。 2 2 4 基于粘性阻尼的动态松弛法 动态松弛法是把非线性静力学问题化为动力学问题求解的一种数值方法式，该方 法的实质是对临界阻尼振动方程进行逐步积分【3 4 1 。为了保证求得准确解，一般都采用 质量阻尼和刚度阻尼来吸收系统的动能，当阻尼系数取值稍小于某一临界值时，系统 的振动将会很快停止，计算趋于收敛。这种带有阻尼项的动态平衡方程，利用有限差 分法按时步求解就是所谓的动态松弛法。该方法的特点是：整个计算过程只需要直接 西南交通大学硕士研究生学位论文第10 页 曼曼曼皇曼曼量曼曼i m _ m m - - 一_ l 曼曼曼 代换，即利用前一时步的函数值计算新的函数值。以单质点的计算为例，离散单元法 的基本运动方程则可修正为 m i i ( t ) + c f 。( t ) + 七。u ( t ) = f ( t )( 2 9 ) 其中，m 是单元的质量；“是位移；f 是时间；c 粘性阻尼系数；砌刚度系数；f 是单 元所受的外荷载。式( 2 9 ) 1 豹i 动态松弛解法就是假定h 缸时刻以前的变量f ( 0 ，甜( 0 ， i i ( t a t ) ，u ( t a t ) 以及u ( t a t ) 等已知，利用中心差分法，可以求得 西( f ) ：u ( t + a t ) - _ 2 u 了( t ) + u ( t 一- a t ) ( 2 - 1 o ) f 矗( f ) ：u ( t + a t ) 石- _ u ( t 一- a t ) ( 2 - 1 1 ) 2 址 将式( 2 1 0 ) 和( 2 1 1 ) 代入式( 2 9 ) ，并整理可得 “o + a t ) = 缸2 厂( f ) + 吒& 一删) 甜。一a t ) + ( 2 m - k a t 2 ) “u ) 1 7 ( 聊+ 丢出) ( 2 - 1 2 ) 从以上可以看出，基于动态松弛法的计算比较简单，也节省计算时间，但是也有 其不足的地方：因为粘性阻尼是人为添加的一个参数，其目的是为了使计算快速收敛， 因而基于动态收敛计算得到的稳定时间并非系统真正的稳定时间。 2 3 几种离散单元法接触模型介绍 2 3 1 弹簧一阻尼器模型 该模型由c u n d a l l 提出阳。如图2 - 4 ( a ) 、图2 4 ( b ) 所示，图中吒为法向刚度， k ，为切向刚度，d 。为法向阻尼，d t 为切向阻尼，以上参数通常根据试验或经验确定， 模型允许颗粒间发生切向滑移及相互分离。 西南交通大学硕士研究生学位论文第11 页 2 3 2 球形颗粒接触模型 矗) 辕向作用b ) 钠自作用 图2 4 颗粒接触模型示意图 除采用弹簧阻尼器给定参数外，也可利用接触力学的成果：法向作用根据h e r t z 理论确定 3 6 】，切向作用根据m i n d l i n e d e r e s i e w i z z 理论确定【3 7 4 2 1 。 假设两个球体发生接触，如图2 - 5 所示。若接触面的半径为口的圆，则接触压力 为 肿，= 等) 2 l 2 3 ， 式中 p ( r 卜接触压应力； 卜接触压应力的合理；口为接触半径。最大接触压应力位于产0 处 3 p p 一2 订 ( 2 - 1 4 ) 一o 西南交通大学硕士研究生学位论文第12 页 图2 5h e r t z 理论模型示意图 接触面变形量与接触面积之间有如下关系 a 2 = r 8 ( 2 - 1 5 ) 式中r 为等效半径，且去= 百1 + 万1 ，r 。、r ：分别为两个球体的半径； 万两颗粒的重叠量。 万= ( x l x 2 ) 2 + ( y l y 2 ) 2 一r l r 2 ( 2 - 1 6 ) 接触面压应力与变形之间的关系为 p ：4 - - - e - er12万32(2-17) 式中e 等效弹性模量，有 丢= 警- 4 卜警 亿一= 二- - ，- i e e 。e 、 式中 置、最颗粒的弹性模量； 卟呸颗粒的泊松比。 如果两个颗粒发生切向接触，首先沿着接触面表面圆周发生相对滑移，其法相与 切向力法相一致，假设两接触表面相对切向位移为口，见图2 - 6 所示， 西南交通大学硕士研究生学位论文第13 页 图2 - 6m i n d l i n e d e r e s i e w i c z 理论模型不葸图 对应的切向力r 可为 t = 8 a g 吼口+ ( 一1 ) 。z ( 1 一o k ) n 式中k = o ，1 ，2 ，分别对应加载，卸载及卸载后重新加载情况； 接触表面的摩擦系数； g 剪切模量当量，且有 g 兰( 2 - v 1 ) g l + ( 2 - 0 2 ) g 2 式中g 1 、g ：颗粒的剪切模量； u 。、d 2 颗粒的泊松比。 吼可以根据式( 2 2 1 ) 确定，如有例 o ；当戈 o ( 图3 - 2 ( c ”时，滑块向左 运动，乃指向右，亩 o 对) ( 3 - 4 ) ，戚+ 缸= 乃 ( 膏 0 ( 3 - 6 ) f x o ) = a 2c o s ( p 。f 一驴2 ) + 之生 戈 o ( 3 7 ) ，【 式中p ：= 一k ，4 、a 2 撕o , 、仍分别取决于各自半周运动中的初始条件。 r h 设t = 0 时，初始位置x = x 0 初始速度量= 童。，且滑块初始位置对应的弹性力大于静 摩擦力，则滑块从右向左开始运动，其位移曲线如图3 3 所示，为了对比，图3 4 给出 具有粘性阻尼振动系统的位移时间曲线。 图3 3c o u l o m b 摩擦阻尼的位移时间曲线 西南交通大学硕士研究生学位论文第19 页 皇鼍i _mm,m i m _ 量曼皇曼! 鼍曼曼曼皇曼曼曼! 寰曼曼 图3 4 粘性阻尼的位移时间曲线 根据以上分析可得，滑块受c o u l o m b 摩擦阻力作用时，在每半个周期内，滑块作 自由振动频率为p 。的简谐振动，但是在每半个周期内速度为零的平衡点都是不同的。 且每经过一个周期的运动，振幅减少阜，按等差级数减少，即x t 曲线的包络线是直 庀 线。当在某一半周期结束的位移所对应的弹性力不足以克服摩擦力时，运动便停止， 最后停止在干阜之间的某一位置，平衡点可根据受力分析求得。相比之下，对粘性 后 阻尼振动中的欠阻尼情形，当f 专0 0 时，欠阻尼系统的位移趋于零，且滑块停在初始位 置。 3 3 单自由度c o u l o m b 摩擦阻尼振动问题的数值解法 3 3 1 计算模型 如图3 5 所示，假定质量为m 的滑块，通过弹簧与墙面相连，弹簧的刚度系数为k ， 接触面的摩擦系数为，所受为外力f 。 图3 5 单个滑块的计算模型 f 一 西南交通大学硕士研究生学位论文第2 0 页 曼曼曼曼笪曼皇曼曼曼曼蔓曼曼曼曼曼曼曼! i 一一一l i elilie i i i i 一一i 。曼曼! 曼曼舅舅! ! ! 曼曼曼曼鼍曼! 皇 在某一时刻，弹簧接触力为 e = 一h( 3 - 8 ) 当速度不为零时，为滑动摩擦力，且 乃= 胪帽( 3 9 ) 当速度为零时，为静摩擦力，为 b = 一k x + f ( 3 1 0 ) 摩擦力的方向与滑块的运动状态、受力有关。如果滑块速度v 0 ，则接触面上的 摩擦力为滑动摩擦力，摩擦力大小根据式( 3 9 ) 确定，方向与速度相反；当速度的绝 对值v = 0 时，滑面上的摩擦力为静摩擦力，且有0 f p m g ，方向与合外力相反。 根据牛顿第二定理，滑块的加速度值可表达为 口：- k x + _ a g n g + f ( 3 - 5 - 1 1 ) 口= 一 i li ， 即式( 2 3 ) 的