刚度系数

刚度系数

乙二醇防冻液配方-海洋污染有哪些

2023年2月19日发(作者：sart)

第ｌ６卷第３期

２００７年９月

计算机辅助工程

ＣＯＭＰＵＴＥＲ ＡＩＤＥＤ ＥＮＧＩＮＥＥＲＩＮＧ

Ｖ０１．１６ Ｎｏ．３

Ｓｅｐｔ．２００７

文章编号：１００６－０８７１（２００７）０３－００６５４３４

颗粒离散元法中阻尼系数、刚度系数和

时步的选取方法

杨 洋， 唐寿高， 王居林

（同济大学航空航天与力学学院，上海２０００９２）

摘要：为解决颗粒离散元法（Ｄｉｓｔｉｎｃｔ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ，ＤＥＭ）中阻尼系数等参数选取困难的问题，

对岩体工程中的ＤＥＭ与颗粒ＤＥＭ进行区分，着重讨论颗粒ＤＥＭ中阻尼系数、刚度系数和时步等

参数的选取方法．这些计算参数的合理选取对保证模拟的真实性具有重要意义．

关键词：颗粒离散元法；阻尼系数；刚度系数；时步

中图分类号：ＴＵ４５２；ＴＢ１１５ 文献标志码：Ａ

Ｓｅｌｅｃｔｉｏｎ ｍｅｔｈｏｄｓ ｏｎ ｄａｍｐｉｎｇ ｃｏｅｆｆｉｃｉｅｎｔ，ｓｔｉｆｆｎｅｓｓ ｃｏｅｆｆｉｃｉｅｎｔ

ａｎｄ ｔｉｍｅ ｓｔｅｐ ｉｎ ｇｒａｎｕｌａｒ ｄｉｓｔｉｎｃｔ ｅｌｅｍｅｎｔ ｍｅｔｈｏｄ

ＹＡＮＧ Ｙａｎｇ，ＴＡＮＧ Ｓｈｏｕｇａｏ，ＷＡＮＧ Ｊｕｌｉｎ

（Ｓｃｈｏｏｌ ｏｆ Ａｅｒｏｓｐａｃｅ Ｅｎｇ．＆Ａｐｐｌｉｅｄ Ｍｅｃｈａｎｉｃｓ，Ｔｏｎ￣ｉ Ｕｎｉｖ．，Ｓｈａｎｇｈａｉ ２０００９２，Ｃｈｉｎａ）

Ａｂｓｔｒａｃｔ：Ｔｏ ｓｏｌｖｅ ｔｈｅ ｐｒｏｂｌｅｍｓ ｏｆ ｓｅｌｅｃｔｉｎｇ ｃｏｍｐｕｔａｔｉｏｎａｌ ｐａｒａｍｅｔｅｒｓ ｉｎ ｇｒａｎｕｌａｒ Ｄｉｓｔｉｎｃｔ Ｅｌｅｍｅｎｔ Ｍｅｔｈ—

ｏｄ（ＤＥＭ），ＤＥＭ ｕｓｅｄ ｉｎ ｒｏｃｋ ｅｎｇｉｎｅｅｒｉｎｇ ａｎｄ ｇｒａｎｕｌａｒ ＤＥＭ ａｒｅ ｄｉｆｆｅｒｅｎｔｉａｔｅｄ．Ｔｈｅ ｓｅｌｅｃｔｉｏｎ ｍｅｔｈｏｄｓ

ｏｎ ｄａｍｐｉｎｇ ｃｏｅｆｆｉｃｉｅｎｔ，ｓｔｉｆｆｎｅｓｓ ｃｏｅｆｆｉｃｉｅｎｔ ａｎｄ ｔｉｍｅ ｓｔｅｐ ａｒｅ ｅｍｐｈａｓｉｚｅｄ．Ｔｈｅ ｒｅａｓｏｎａｂｌｅ ｓｅｌｅｃｔｉｏｎ ｏｆ

ｔｈｅｓｅ ｃｏｍｐｕｔａｔｉｏｎａｌ ｐａｒａｍｅｔｅｒｓ ｐｌａｙｓ ａｎ ｉｍｐｏｒｔａｎｔ ｒｏｌｅ ｉｎ ｇｕａｒａｎｔｅｅｉｎｇ ａｕｔｈｅｎｔｉｃｉｔｙ ｏｆ ｓｉｍｕｌａｔｉｏｎ．

Ｋｅｙ ｗｏｒｄｓ：ｇｒａｎｕｌａｒ ｄｉｓｔｉｎｃｔ ｅｌｅｍｅｎｔ ｍｅｔｈｏｄ；ｄａｍｐｉｎｇ ｃｏｅｆｆｉｃｉｅｎｔ；ｓｔｉｆｆｎｅｓｓ ｃｏｅｆｆｉｃｉｅｎｔ；ｔｉｍｅ ｓｔｅｐ

０ 引 言

离散元法（Ｄｉｓｔｉｎｃｔ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ，ＤＥＭ）是由

ＣＵＮＤＡＬＬ 提出的１种处理非连续介质问题的数

值模拟方法，其理论基础是结合不同本构关系的牛

顿第二定律，采用动态松弛法求解方程．

ＤＥＭ自问世以来，其主要应用领域集中在岩体

工程和粉体（颗粒散体）工程．首先，在岩体计算力

学方面，由于离散单元能更真实地表达节理岩体的

几何特点，便于处理所有非线性变形和破坏都集中

在节理面上的岩体破坏问题，被广泛应用于模拟边

坡、滑坡和节理岩体地下水渗流等力学过程．其次，

在粉体工程方面，颗粒离散元被广泛应用于粉体在

复杂物理场作用下的复杂动力学行为的研究和多相

混合材料介质或具有复杂结构的材料力学特性研究

中．它涉及到粉末加工、研磨技术、混合搅拌等工业

加工领域以及粮食等颗粒离散体的仓储和运输等生

产实际领域。

岩体工程中的ＤＥＭ与颗粒ＤＥＭ并无本质不

同，但在接触处理以及一些概念的认识上有一定区

别。例如，在节理岩体问题中，单元之间总是处于相

互接触或存在接触一断开的过程，均可视为准静态

收稿日期：２００７－０３－１２修回日期：２００７－０７－１６

作者简介：杨洋（１９８１－），男，湖南安化人，硕士，研究方向为离散元法的数值模拟，（Ｅ－ｍａｉｌ）ｙｙａｎｇ１６７＠ｇｍａｉｌ．ｔｏｍ；

唐寿高（１９５２－），男，上海人，教授，博士，研究方向为计算力学，（Ｅ－ｍａｉｌ）ｔａｎｇｓｇ＠ｍａｒｌ、ｔｏｎｇｊｉ、ｅｄｕ．ｃｎ

维普资讯

计算机辅助工程 ２００７皇

情况，在此基础上引人动态松弛法 将该准静态问

题化为动力学问题进行求解．动态松弛法要求选取

合适的阻尼，使函数收敛于静态值．在颗粒体问题

中，颗粒间并不一定总存在接触，颗粒体间的相互碰

撞也表现为动态的过程，此时采用动态松弛法进行

求解并非为了得到静态值，而是为了引人阻尼系数

以提供耗能装置，达到最大程度的模拟效果．

本文旨在对颗粒ＤＥＭ中阻尼等计算参数的选

取方法进行阐述，有关ＤＥＭ原理的详细论述可参考

文献［３］．

１ 阻尼系数选取

颗粒ＤＥＭ中阻尼系数的选取可参考连续介质

中阻尼的取法，引人工程中的黏性阻尼概念，采用

Ｒａｙｌｅｉｇｈ线性比例阻尼．Ｒａｙｌｅｉｇｈ线性比例阻尼可以

表示为

［Ｃ］＝ ］＋ Ｋ］ （１）

式中：［ｃ］为阻尼矩阵；［ ］为质量阻尼矩阵；［ ］

为刚度阻尼矩阵； 和 分别为质量阻尼比例系数

和刚度阻尼比例系数，可由下式确定

＝ ｉ ｃｃＪ （２）

＝ ／ｔｏ （３）

式中： ｉ 和∞ ｊ 分别为系统的最小临界阻尼系数和

最小圆频率．对于弹性连续介质，可按下式求得

２１Ｔ厂— 丁 一 …

（－Ｏｍｉｎ √ ：

常用的系统振动阻尼比 的确定方法有半功率

法和对数减量法等．

如前所述，Ｒａｙｌｅｉｇｈ阻尼理论适用于连续介质

系统，不完全适用于颗粒体这样的非连续介质系统，

因为非连续介质系统随着单元之间的滑移或分离，

其振型不确定，但阻尼却仍然存在，并可以用图ｌ所

示的物理模型解释．可以想象图中质量阻尼ｄ 为把

整个系统浸泡在黏性液体中，在物理意义上等价于

用黏性活塞将颗粒单元与一不动点相连，使块体单

元的绝对运动受到阻尼．刚度阻尼ｄ 在物理意义上

等价于用黏性活塞把两个接触块体相连，使颗粒单

元之间的相对运动受到阻尼．

当颗粒之间接触完全脱离，即不存在颗粒之间

的相互接触时，阻尼不再存在，或者将此时的阻尼理

解为颗粒在空气中受到的质量阻尼．所以，在颗粒

ＤＥＭ中，实际存在一个变阻尼的概念，包含至少两

套阻尼，即接触时的质量阻尼加刚度阻尼和无接触

时的空气质量阻尼．

对于连续介质来说，其振型、最小圆频率 和

（ｂ）

图１阻尼的物理模型

最小临界阻尼系数 ｉ 等能够经过计算与实验得

到．但是，对于非连续介质，由于其振型不确定，只能

用试算的办法确定这些参数进而计算阻尼系数．颗

粒ＤＥＭ中引人阻尼系数是为了提供耗能装置，并非

为了得到准静态解，因此，阻尼系数的选取具有一定

的灵活性，以满足最大程度模拟为原则．

２ 刚度系数选取

对刚度系数的考虑见图２，颗粒体Ａ与颗粒体

Ｂ存在两个角边接触，接触力分别为Ｆ。和Ｆ２，对于

块体Ａ有平衡方程

Ｆｌ＋Ｆ２＝Ｐ ＝ｐｌｈｇ （５）

ｒ．．— —１

图２ ｋ 计算简图

式中：Ｖ和ｇ分别为颗粒体Ａ的体积和重力加速度，

块体厚度为 设两个角边接触的稳定叠合量为△ ，

则式（５）可写为

２，４ ｋ ＝ｐｌｈｇ （６）

设颗粒体的弹性模量为Ｅ，可导出颗粒体Ａ顶

部的位移量

Ａ＝ ｄ＾＝ （７）

令Ａ ＝Ａ，则有

一 ＝三 一ｆ

维普资讯

第３期 杨 洋，等：颗粒离散元法中阻尼系数、刚度系数和时步的选取方法 ６７

ｋ ＝÷Ｅ （８）

，

对于ｋ ，采用下式计算：

ｋ ＝ｋ．／［２（１＋ ）］ （９）

式中， 为单元材料的泊松比．

模拟实践表明， 和ｋ 取值过小时，颗粒之间

的嵌入量将过大而引起计算失真；另一方面，ｋ 和

ｋ 取值过大时，则由于弹簧过刚而引起单元的反弹

失稳．计算收敛时，刚度系数ｋ 和ｋ 大到一定程度

后，对计算结果影响不大．

３ 时步选取

时步计算的理论基础是求解单自由度有阻尼弹

性体系的中心差分格式下的Ｉ临界时步△￡．对于动力

方程

ｍｉｉ＋ｋｕ＝０ （１０）

式中：ｍ为颗粒单元质量；Ｕ为位移；ｋ为刚度系数．

由中心差分法，ｉｉ可写成：

ｉｉ＝［Ｕ（ｔ＋△￡）一２ （ｔ）＋Ｕ（ｔ—Ａｔ）］／（△￡）

将上式代入式（１０）中，有

ｍ［Ｕ（ｔ＋△￡）一２ （ｔ）＋Ｕ（ｔ—Ａｔ）］／

（△￡） ＋ （ｔ）＝０

或者

Ｕ（ｔ＋△￡）＋［ｋ（Ａｔ）。／ｍ一２］Ｕ（ｔ）＋

Ｕ（ｔ—Ａｔ）＝０ （１１）

根据差分理论，方程（１１）的解为

（￡）：『２一－－ｋ（ａｆ） ±

√（ ） ｃ△ 一４ ｃ△ 。］／２

为了使解具有振荡性，Ｕ（ｔ）必须是复数，即要求：

（ ） （ －－Ｌｍ（ａ ＜０

或

△ ＜２√詈 ２ （１２）

又由∞ ＝２７ｒ／Ｔ（Ｔ是固有振动周期），可把上式写成：

△￡＜一Ｔ （１３）

１Ｔ

理论证明，系统的最小固有振动周期总是大于

其中任何一个单元的最小固有振动周期 ｉ ，将后

者用于时步计算，其结果是安全的．因此，在ＤＥＭ计

算中通常取时步

Ａｔ＜Ｔｍｉ ／１０ （１４）

式中

ｎ ｍｉｎ（＾√／等） （１５）

式中，ｍｉｎ（）表示取最小值．

由推导可知，采用上述方法计算的时步能够达

到足够小，可以保证颗粒之间的接触过程得到充分

模拟，不会出现这个时步颗粒之间刚刚开始接触，下

个时步颗粒间的接触就反弹开了的现象，保证了接

触模拟的真实性．

４ 算 例

下面给出采用本文作者编制的颗粒ＤＥＭ筒仓

计算程序ＳＩＳＯＬＶ．２¨４］，对某大型筒仓的装、卸料过

程进行模拟的算例．对原６０ ｍ直径、２０ ｍ仓高的筒

仓按２５：３缩小建立模型，模型尺寸见图３．模拟中

采用的计算参数见表１．

｝＿ ｍ 卜 卜 一十４ｓ｝＿ｔｍ

表１模拟中采用的计算参数

图３简仓模型几何尺寸

模型初始状态见图４，物料在重力的作用下压

紧密实，经过１．０×１０ 个时步的运算，达到静止．此

时模型状态见图５．

物料达到静止后，开始开仓卸料．图６—８描绘

开仓１５ ０００个时步后模型的运动状态、速度场和压

力场．

在图７中，每根线段表示该点处颗粒的速度，线

段长度和方向分别表示速度的大小和方向．在图８

中，每根线段表示该点处的作用力，线段长度和方向

分别表示压力的大小和方向．通过对模型作速度场

和压力场图，能够很直观地看到整个系统的速度和

压力分布情况．

Ｔ ● ●上

维普资讯

计算机辅助工程 ２００７车

图４筒仓模型初始状态

图６开仓１５ ０００个时步后模型运动状态

图７开仓１５ ０００个时间步后模型速度场

图８开仓１５ ０００个时间步后模型压力场

５ 讨论

颗粒ＤＥＭ看似简单，其实却很难．如何选取上

述几个参数对于初学者是很棘手的问题．要得到正

确的模拟结果，需要在深入理解某些相关概念的基

础上通过试算得到阻尼等计算参数，只有选取合理

的计算参数才能保证模拟的真实性．

参考文献：

［１］ＣＵＮＤＡＬＬ Ｐ Ａ．Ａ ｃｏｍｐｕｔｅｒ ｍｏｄｅｌ ｆｏｒ ｓｉｍｕｌａｔｉｎｇ ｐｒｏｇｒｅｓｓｉｖｅ ｌａｒｇｅ ｓｃａｌｅ ｍｏｖｅｍｅｎｔｓ ｉｎ ｂｌｏｃｋｙ ｒｏｃｋ ｓｙｓｔｅｍ［Ｃ］／／Ｐｒｏｃ Ｓｙｍｐｏｓｉｕｍ ｏｎ Ｉｎｔｅｒｎａｔｉｏｎａｌ

Ｓｏｃｉｅｔｙ ｏｆ Ｒｏｃｋ Ｍｅｃｈａｎｉｃｓ，Ｆｒａｎｃｅ：Ｎａｎｃｙ，１９７１（１）：８－１２．

［２］ＯＴＹＥＲ Ｊ Ｒ Ｈ，ＣＡＳＳＥＬＬ Ａ Ｃ，ＨＯＢＢＳ Ｒ Ｅ．Ｄｙｎａｍｉｃ ｒｅｌａｘａｔｉｏｎ［Ｃ］／／Ｉｎｓｔｉｔｕｔｉｏｎ ｏｆ Ｃｉｖｉｌ Ｅｎ ｎｅｅｒｓ Ｐｒｏｃ，１９６６，３５（４）：６３３４５５６．

［３］ 王泳嘉，邢纪波．离散单元法及其在岩－＋－２学中的应用［Ｍ］．沈阳：东北工学院出版社，１９９１．

［４］ 杨洋．水泥熟料多点卸料的离散元法模拟［Ｄ］．上海：同济大学，２００７．

［５］ 汪泽森，董玉德．自适应参数化ＣＡＤ图块数据库系统的设计［Ｊ］．计算机辅助工程，２００５，１４（２）：１４．

（编辑廖粤新）

２００７产品现代设计与技术创新论坛举行

２００７产品现代设计与技术创新论坛于２００７年９月１Ｏ一１ １日在上海交通大

学举行，国内外著名企业和中小企业的代表、高等院校的代表共一百多人参加了

本届论坛．本界论坛的宗旨是推动现代设计理论、方法与工具在中国企业的应

用，增强中国企业的创新能力，主题为”创新模式设计方法知识服务评价体

系”．论坛上中国工程院谢友柏和郭孔辉院士、美国南加州大学Ｓｔｅｐｈｅｎ Ｃ—Ｙ Ｌｕ

教授等知名专家学者和企业代表进行了主题报告，之后代表们进行了热烈的讨

论．

维普资讯