离散元课件(二forcopy)分解

二基本原理-球形颗粒元离散元法岩土工程研究所刘军根据离散化模型中所采用的单元种类分别介绍离散元法的基本原理：颗粒元•二维圆盘单元•三维圆球单元块体元•多边形单元•多面体单元二基本原理基本假设假定速度和加速度在每个时间步长内为常量；选取的时间步长应该足够小以至于在单个时间步长内扰动的传播不会超过当前与之相邻的粒子。二基本原理-球形颗粒元离散元法运动描述处于一个理想散体中的任意一个颗粒,具有6个自由度,3个平动自由度与三个转动自由度,可通过Newton第二定律分别描述。二基本原理-球形颗粒元离散元法运动描述平动方程：二基本原理-球形颗粒元离散元法gFFVikjijdijciimtmi1,,dd式中，与分别为颗粒的质量和速度。为时间，为颗粒的重力，与分别为颗粒与的接触力与粘性接触阻尼力，为所有与颗粒接触的颗粒总数。imiVitgimijc,Fijd,Fijik运动描述接触力的分解：颗粒与间的接触力可分解为法向与切向接触力，即二基本原理-球形颗粒元离散元法jiijctijcnijc,,,FFF同理，粘性接触阻尼力也可分解为法向与切向分量形式，即ijdtijdnijd,,,FFF运动描述接触力产生的力矩：二基本原理-球形颗粒元离散元法颗粒间的接触力作用在两个颗粒的接触点上，而不是作用在颗粒的中心，所以这些接触力（除法向接触力外）将会对颗粒产生力矩，ijcn,FiTijdtijctii,,FFRT式中，为从颗粒的质心指向接触点的矢量，其幅值为（颗粒的半径）。iRiiR运动描述转动方程：二基本原理-球形颗粒元离散元法转动方程可以表示为ikjiiitI1ddTω式中，与分别为颗粒的转动惯量与角速度，对于球形颗粒为iIiωi252iiiRmIiI接触模型综述：二基本原理-球形颗粒元离散元法关于接触力的计算模型已有大量的研究成果，目前仍旧是一个活跃的研究领域，特别是对于切向力的计算方法。对于理想散体颗粒（无粘连）：采用Hertz理论描述法向作用，而采用Mindlin与Deresiewicz理论描述切向作用；对于存在粘连的散体颗粒：法向接触力根据在Hertz理论基础上考虑粘连力的JKR（Johnson-Kendall-Roberts）理论确定，切向接触力增量则根据把Savkoor和Briggs理论与Mindlin和Deresiewicz理论相结合形成的理论确定。接触模型两个处于接触颗粒单位法向和切向向量：二基本原理-球形颗粒元离散元法iiR/Rn单位法向向量单位切向向量nnnntijijijijVVVV单位切向量之所以通过两个颗粒的相对速度来计算，是因为接触力与粘性阻尼力的方向与相对速度的方向相同。接触模型两个处于接触颗粒接触点的相对速度：二基本原理-球形颗粒元离散元法法向相对速度为iijjRωRωVVVijijnnijijn,VV切向相对速度为nnijijVVVijt,nnijijt,VV或者写为接触模型法向接触力计算模型—Hertz模型：二基本原理-球形颗粒元离散元法为颗粒i与j接触时的侵入深度n2/3**,34nijcnREF式中)1(22*vEEjiRRR11*nδijjinRRRR接触模型法向接触力计算模型—Cundall模型：二基本原理-球形颗粒元离散元法nnnijcnk,F式中，为法向弹簧刚度。nk接触模型法向接触力计算模型—法向粘性接触阻尼力：二基本原理-球形颗粒元离散元法式中，为法向粘性接触阻尼系数。nnijVFnijdnc,nc接触模型切向接触力计算模型—综述：二基本原理-球形颗粒元离散元法处于接触中的两个颗粒的切向作用，从本质上讲，是一种摩擦行为，按照摩擦机理，摩擦力包括：滑动摩擦、滚动摩擦与静摩擦，其中滑动摩擦与静摩擦属于切向摩擦力；滚动摩擦是由于法向接触应力的不均匀分布产生的。介绍两个切向接触力模型：Coulomb准则Mindlin与Deresiewicz切向接触力模型接触模型切向接触力计算模型—Coulomb准则：二基本原理-球形颗粒元离散元法在离散元模拟中，一般用Coulomb准则这种简单的形式描述，静摩擦的详细刻画需要涉及切向位移甚至可能要考虑时间依赖效应。式中，为静摩擦系数，切向摩擦力的方向为与相对滑动的趋势相反。sμijcnsijctijcnsijcnsijctijctijct,,,,,,,,,FFFFFFF接触模型切向接触力模型—Mindlin与Deresiewicz模型：二基本原理-球形颗粒元离散元法tmaxt,maxt,ttijcnsijctδδμδδδFF23,,,min11式中为颗粒与间的累积切向位移矢量tδijnsmaxt,δvvμδ122tNtNtnnijVδδ1接触模型切向接触力模型—阻尼力：二基本原理-球形颗粒元离散元法式中，为切向粘性接触阻尼系数。nnijVFtijdtc,tc二基本原理-球形颗粒元离散元法计算模型总结运动方程接触力的计算•法向接触力tmaxt,maxt,ttijcnsijctδδμδδδFF23,,,min11gFFVikjijdijciimtmi1,,ddikjiiitI1ddTωn2/3**,34nijcnREFnnnijcnk,Fijcnsijctijcnsijcnsijctijctijct,,,,,,,,,FFFFFFFnnijVFnijdnc,nnijVFtijdtc,•切向接触力二基本原理-球形颗粒元离散元法运动方程的求解gFFVikjijdijciimtmi1,,ddikjiiitI1ddTω一般采用两种方法求解运动方程：中心差分法Verlet积分法hXXXnnn2121hXXXnnn211二基本原理-球形颗粒元离散元法运动方程的求解中心差分法运动方程可由Verlet显式积分求解。通过积分可获得粒子的新位置，积分时需要粒子的当前及上一步长的位置数据，而不需要粒子的速度数据。二基本原理-球形颗粒元离散元法运动方程的求解Verlet积分法hXXXnnn211hXXXnnn21122242hXXXnnn2112hXXXnnn二基本原理-球形颗粒元离散元法运动方程的求解Verlet积分法定义：从而导出Verlet方程为2112hXXXXnnnn接触发现算法在一个由众多颗粒组成的体系中，直接判别颗粒是否接触需要耗费大量的计算时间，因而，为了节约计算时间，提高计算效率，一般不直接判别任意两个颗粒间是否存在接触，而是分两个步骤判别颗粒间的接触是否存在：首先，对一个颗粒，判别其潜在的邻居个数，然后，准确确定该颗粒与每个邻居是否接触。虽然在确定邻居数目时也要耗费一定的计算时间，但是仍旧比逐个准确判别颗粒间接触是否存在要节约时间。因而，接触发现算法的效率在多颗粒体系力学行为模拟中至关重要。二基本原理-球形颗粒元离散元法接触发现算法介绍三种针对球形颗粒的接触发现算法：Verlet邻居目录法连接单元法边界盒法二基本原理-球形颗粒元离散元法接触发现算法Verlet邻居目录法二基本原理-球形颗粒元离散元法当需要判别体系中某个颗粒的邻居数量时，在该粒子周围构建一个球（称之为参考球，称该颗粒为核心颗粒），参考球半径为体系中最大粒子半径的若干倍，那么参考球所包围的所有粒子为该球中心粒子的邻居。参考球半径的选取取决于粒子的运动速度及体系中粒子的密度。对于每个粒子，都可生成一个邻居粒子的目录。为了得到邻居目录，对每一个粒子而言，所有标号大于该粒子的粒子都必须被检验，判断是否位于该粒子的参考球中，而对于标号小于该粒子标号的粒子则没有必要被检验，因为邻居是互相的，没有必要对一个邻居对检验两次。Verlet邻居目录法及粒子存储目录二基本原理-球形颗粒元离散元法接触发现算法Verlet邻居目录法二基本原理-球形颗粒元离散元法接触发现算法Verlet邻居目录法对于n个颗粒组成的体系，用Verlet邻居目录法需要次计算，也就是说计算次数仍旧为量级。然而，并不需要在每个时间步长上都对邻居目录进行更新。更新的频率取决于体系中粒子的密度、粒子的运动速度以及参考球的尺寸。参考球的半径也可以根据颗粒体系的稠密程度及运动速度进行动态调整，并且参考球半径与邻居目录的更新频率呈反比关系，参考球半径越小，邻居目录的更新频率越高；但是参考球半径越大，则有更多的粒子位于球体内，所以判别是否为邻居就需要较长的时间。2/1nn2no接触发现算法连接单元法二基本原理-球形颗粒元离散元法将颗粒体系所占据的空间划分成规则的网格，对于三维问题，可以划分为个立方体单元，对于二维体系，则划分为个正方形单元，对于颗粒体系所占据空间形状不规则时，也可采用其他形状单元划分。但是所有单元的尺寸必须大于粒子的尺寸。与Verlet邻居目录法的主要区别在于：相邻单元法中的单元不依附于粒子，单元不随粒子的运动而运动。如果粒子就当前的位置被分配到某个单元，显然，只有在同一个单元或直接相邻单元内的粒子间才可能发生相互作用，也就是说，只有相邻单元内的粒子才能成为邻居。mmmmm接触发现算法连接单元法二基本原理-球形颗粒元离散元法例如，对一个二维体系而言，只有在9个不同的单元内可能包含邻居粒子，对于三维体系，则只有在27个不同单元内可能包含邻居粒子。与前面介绍的Verlet邻居目录法相同，每个粒子对只需要检验一次，这样，没必要对9个单元都进行检验，只需检验中心单元及邻居单元的一半即可，即在二维情况下，只需检验5个单元，在三维情况下，只需检验14个单元。接触发现算法连接单元法二基本原理-球形颗粒元离散元法2D体系中的连接单元法接触发现算法边界盒法二基本原理-球形颗粒元离散元法这个方法与前述两种方法不同。首先，在每一个粒子周围构建一个边界盒，边界盒的尺寸按这样的方式选取：使每个粒子刚好放进它的边界盒内。边界盒的边为直线，并且与体系的坐标轴平行。在判别颗粒的邻居时，把边界盒投影到体系的坐标轴上。通过边界盒在坐标轴上投影的起点和终点来判别是否为邻居。接触发现算法边界盒法二基本原理-球形颗粒元离散元法a围绕在每个粒子周围的边界盒b两个不同时刻粒子的边界盒在轴上的投影在判别颗粒的邻居时，把边界盒投影到体系的坐标轴上。例如，图a表明了粒子及边界盒的位置，图b为图a中的边界盒在体系x轴的投影。通过边界盒在坐标轴上投影的起点和终点来判别是否为邻居，出于这个原因，投影的起点和终点序列被存储在目录中。接触发现算法边界盒法二基本原理-球形颗粒元离散元法对于3D体系，必须把边界盒子在三个坐标轴上投影，所以生成三个目录。如果一个粒子的边界盒在某个坐标轴上投影的起点和终点间包含另外一个粒子边界盒投影的起点、终点或起点和终点，就说明这两个粒子的边界盒在该坐标轴上的投影发生了重叠。如果两个边界盒的投影在每一个坐标轴上都发生重叠，那么就说明这两个边界盒发生了接触。接触发现算法边界盒法二基本原理-球形颗粒元离散元法存储目录的更新检查一个边界盒在每一个坐标轴的投影是否在另外一个边界盒投影的起点和终点之间仍旧需要花费很多计算时间。但是，尽管这些目录在每个计算时间步长上都必须被更新，所用的计算时间仍旧可以被减少到与体系中粒子数量成比例的量级，因为在每一个新步长上所需要做的只是对旧目录的更新，并且对旧目录的更新只是对旧目录的重新分类，因而，更新的过程非常简单，只需要对旧目录进行顺序地检查，判断是否在顺序上有新的变化即可。这种变化仅仅是位置的变化。二基本原理-球形颗粒元离散元法算法设计二基本原理-球形颗粒元离散元法计算所需参数只有当选取的时间步长小于临界时间步长时，数值计算才是稳定的，这是模型采用显式积分的