一种基于距离势函数三维任意凸多边形块体离散单元法制造技术

本发明专利技术公开了一种基于距离势函数三维任意凸多边形块体离散单元法，将研究对象离散为块体单元体系；根据本发明专利技术提供的非均匀块体离散单元接触检测方法，对所有块体单元进行接触检测；定义距离势函数，对相互接触的块体单元进行法向接触力计算，得到作用于块体单元的法向接触力及其引起的力矩；基于本发明专利技术提供的切向接触力方向矢量计算方法，根据切向相对位移的增量，计算得到作用于块体单元的切向接触力及其引起的力矩；由刚体运动方程求解块体单元的加速度与角加速度；采用velocity verlet算法计算块体单元的速度与位移，根据各块体单元的运动和相对位置，描述整个介质的变形和运动状态。

【技术实现步骤摘要】

本专利技术涉及一种块体离散单元法，具体涉及一种基于距离势函数三维任意凸多边形块体离散单元法，属于块体离散元

技术介绍
离散单元法的思想源于较早的分子动力学，1971年由美国的CundallP.A.教授首先提出用以模拟非连续介质力学行为的数值分析方法。这种方法是将研究对象划分为大量的离散块体(颗粒)，每个块体拥有相应的自由度。块体之间通过接触力、力矩相互作用，同时施加合适的力或位移的边界条件。在给定初始条件及外力作用下，基于牛顿第二运动定律并采用显示时步步进的动态松弛法，确定下一时刻块体的运动状态。离散单元法原理直观简单，以每个块体的运动为基础，构建整个系统的运动情况。能够模拟研究物体在动态条件下的变形和破坏过程，特别是大变形以及断裂等。目前离散单元法已经在散体力学、岩土工程和地质工程等领域得到广泛应用并取得重要成果。目前，离散单元法主要分为两种。一是由Cundall教授提出传统离散单元法。该方法在计算中可能存在能量不守恒的问题；接触检测判断需考虑角-角接触，角-面接触，面-面接触等不同的块体接触形式，过程繁琐；无法处理角-角接触的情况，需采用圆角化进行规避，过程复杂；接触检测方面，目前普遍采用的方法是公共面法，但是公共面法仍然存在很多问题，例如公共面的正确选取、唯一性以及旋转时的迭代误差等。二是由英国A.MUNJIZA教授提出有限离散单元法，通过将研究对象划分为大小均匀的四面体块体单元，并以单元形心为基础建立势函数定义以此计算单元间的接触力；接触检测方面提出了针对大小均匀单元接触检测的线性检测方法—NBS检测方法。与Cundall提出的传统离散单元法相比，有限离散单元法很好地规避了不同接触形式的接触检测问题，解决了传统方法不能处理角-角接触的问题，且整个过程满足能量守恒的要求。但仍存在一些问题：应用大小均匀的四面体单元，一方面模型与实际情况不符，另一方面在实际应用时，均一化的单元尺寸以及最简单的单元形式会大大增加划分块体单元的数量，降低计算效率；单元的刚度以及法向接触力的计算受单元形式影响；接触力计算中没有考虑切向接触力；由于任意形状的凸多面体单元的形心无法将单元分割为体积相等的子单元，基于单元形心的势函数定义无法解决任意形状的凸多面体块体接触力的情况。
技术实现思路
本专利技术所要解决的技术问题是提供一种基于距离势函数三维任意凸多边形块体离散单元法，解决三维任意凸多面体块体离散单元接触力计算的问题，改善有限离散单元法的不足。本专利技术的方法解决了现有技术中单元尺寸均一、形式单一、单元刚度及接触力的计算受单元形式影响、没有考虑切向接触力以及无法应用于任意多面体单元等技术问题。本专利技术为解决上述技术问题采用以下技术方案：本专利技术提供了一种基于距离势函数三维任意凸多边形块体离散单元法，包括以下步骤：步骤一，将研究对象离散为多面体块体单元体系，并根据研究对象的大小，确定计算区域；步骤二，确定计算时间步长；步骤三，在当前时间步，基于非均匀块体离散单元接触检测法对所有块体单元进行接触检测，并依据块体单元形心之间的距离、块体单元各顶点与其对应的形心之间的最大距离确定每个块体单元及与之接触的块体单元；步骤四，根据步骤三的检测结果，对相互接触的块体单元进行接触力计算，当循环到某个块体单元时，以此块体单元为目标块体单元，与此块体单元相接触的块体单元为接触块体单元，进一步进行接触判断，并基于块体离散单元距离势函数的定义，计算当前时间步作用于目标块体的法向接触力、切向接触力以及力矩；步骤五，重复步骤四计算当前时间步内每个块体单元所受的法向接触力、切向接触力和力矩；步骤六，根据刚体运动方程，计算出每个块体单元当前时间步的加速度，再由velocityverlet算法计算出每个块体单元下一时间步的速度和位移；步骤七，根据步骤六中块体单元的位移，更新每个块体单元顶点、形心的坐标，完成当前时间步的计算；步骤八，重复步骤三至七计算下一时间步，直至计算完所有时间步。作为本专利技术的进一步优化方案，步骤二中的计算时间步长须满足为：其中，Δt是计算时间步长；ξ是系统的阻尼比，m是块体单元质量，c是阻尼系数，k是刚度系数。作为本专利技术的进一步优化方案，步骤三中，基于非均匀块体离散单元接触检测法对所有块体单元进行接触检测的具体步骤包括：步骤3.1，对每个块体单元，计算其最大外接球半径；步骤3.2，所有块体单元根据其最大外接球半径，按从大到小的顺序分成n组，具体为：1)第1组中块体单元的最大外接球半径为d1，d1的取值范围为：D≤d1＜D/α；2)第2组中块体单元的最大外接球半径为d2，d2的取值范围为：D/α≤d2＜D/α2；3)第3组中块体单元的最大外接球半径为d3，d3的取值范围为：D/α2≤d3＜D/α3；4)以此类推，第n-1组中块体单元的最大外接球半径为dn-1，dn-1的取值范围为：D/αn-2≤dn-1＜D/αn-1；5)第n组中块体单元的最大外接球半径为dn，则dn的取值范围为：dn≥D/αn-1；其中，n＝int[log(D/d)/logα+0.416]+1,D为所有块体单元中最大块体单元外接球半径，d为所有块体单元中最小块体单元外接球半径，α＝2；步骤3.3，对第1组的块体单元与所有块体单元进行接触检测，具体步骤包括：1)取lmax＝2d1，将计算区域划分为以lmax为边长的正方体网格；2)根据公式：将所有块体单元映射到网格上，其中，(xK,yK,zK)为第K个块体单元的形心坐标，K＝1,2,…,N，N为块体单元个数，xmin、ymin、zmin分别为计算区域在x、y、z方向的最小值；3)对每个块体单元与其周边所有相邻网格内块体单元进行接触检测，当循环到某个块体单元时，即以此块体单元所在网格为中心，判断此块体单元是否与该网格内部以及相邻网格内部的块体单元接触，其中，判断此块体单元是否与该网格内部以及相邻网格内部的块体单元接触的具体过程为：3-1)分别计算两个块体单元的形心之间的距离dist、两个块体单元各顶点与其对应的形心之间的最大距离，并计算两个块体单元各顶点与其对应的形心之间的最大距离的平均值lavmax；3-2)比较dist与lavmax，若dist≤2lavmax，则判断此两个块体单元接触，否则判断此两个块体单元不接触；步骤3.4，将第一组的块体单元从所有块体单元之中去除，对第2组块体单元与剩余第2组到第n组的块体单元进行接触检测，具体步骤包括：1)取lmax＝2d2，将计算区域划分为以lmax为边长的正方体网格；2)根据步骤3.3中2)至3)的方法，将第2组到第n组的块体单元映射到网格上，并对第2组块体单元与剩余第2组到第n组的块体单元进行接触检测；步骤3.5，以此类推，将已进行接触检测的块体单元组从所有块体单元中去除，并进行下一组块体单元与剩余所有块体单元组的接触检测，直至完成n组块体的接触检测。作为本专利技术的进一步优化方案，步骤四中，目标块体单元的法向接触力、切向接触力以及力矩具体计算步骤如下：步骤4.1，确定目标块体单元的最大内切球直径；步骤4.2，定义距离势函数，循环计算当前时间步每个与目标块体单元相接触的块体单元作用于目标块体的接触力和力矩；步骤4.3，计算有与目标块体单元相接触的块体单本文档来自技高网...

【技术保护点】
一种基于距离势函数三维任意凸多边形块体离散单元法，其特征在于，包括以下步骤：步骤一，将研究对象离散为多面体块体单元体系，并根据研究对象的大小，确定计算区域；步骤二，确定计算时间步长；步骤三，在当前时间步，基于非均匀块体离散单元接触检测法对所有块体单元进行接触检测，并依据块体单元形心之间的距离、块体单元各顶点与其对应的形心之间的最大距离确定每个块体单元及与之接触的块体单元；步骤四，根据步骤三的检测结果，对相互接触的块体单元进行接触力计算，当循环到某个块体单元时，以此块体单元为目标块体单元，与此块体单元相接触的块体单元为接触块体单元，进一步进行接触判断，并基于块体离散单元距离势函数的定义，计算当前时间步作用于目标块体的法向接触力、切向接触力以及力矩；步骤五，重复步骤四计算当前时间步内每个块体单元所受的法向接触力、切向接触力和力矩；步骤六，根据刚体运动方程，计算出每个块体单元当前时间步的加速度，再由velocityverlet算法计算出每个块体单元下一时间步的速度和位移；步骤七，根据步骤六中块体单元的位移，更新每个块体单元顶点、形心的坐标，完成当前时间步的计算；步骤八，重复步骤三至七计算下一时间步，直至计算完所有时间步。...

【技术特征摘要】
1.一种基于距离势函数三维任意凸多边形块体离散单元法，其特征在于，包括以下步骤：步骤一，将研究对象离散为多面体块体单元体系，并根据研究对象的大小，确定计算区域；步骤二，确定计算时间步长；步骤三，在当前时间步，基于非均匀块体离散单元接触检测法对所有块体单元进行接触检测，并依据块体单元形心之间的距离、块体单元各顶点与其对应的形心之间的最大距离确定每个块体单元及与之接触的块体单元；步骤四，根据步骤三的检测结果，对相互接触的块体单元进行接触力计算，当循环到某个块体单元时，以此块体单元为目标块体单元，与此块体单元相接触的块体单元为接触块体单元，进一步进行接触判断，并基于块体离散单元距离势函数的定义，计算当前时间步作用于目标块体的法向接触力、切向接触力以及力矩；步骤五，重复步骤四计算当前时间步内每个块体单元所受的法向接触力、切向接触力和力矩；步骤六，根据刚体运动方程，计算出每个块体单元当前时间步的加速度，再由velocityverlet算法计算出每个块体单元下一时间步的速度和位移；步骤七，根据步骤六中块体单元的位移，更新每个块体单元顶点、形心的坐标，完成当前时间步的计算；步骤八，重复步骤三至七计算下一时间步，直至计算完所有时间步。2.根据权利要求1所述的一种基于距离势函数三维任意凸多边形块体离散单元法，其特征在于，步骤二中的计算时间步长须满足为：Δt≤2·mk·(1+ξ2-ξ)]]>其中，Δt是计算时间步长；ξ是系统的阻尼比，m是块体单元质量，c是阻尼系数，k是刚度系数。3.根据权利要求1所述的一种基于距离势函数三维任意凸多边形块体离散单元法，其特征在于，步骤三中，基于非均匀块体离散单元接触检测法对所有块体单元进行接触检测的具体步骤包括：步骤3.1，对每个块体单元，计算其最大外接球半径；步骤3.2，所有块体单元根据其最大外接球半径，按从大到小的顺序分成n组，具体为：1)第1组中块体单元的最大外接球半径为d1，d1的取值范围为：D≤d1＜D/α；2)第2组中块体单元的最大外接球半径为d2，d2的取值范围为：D/α≤d2＜D/α2；3)第3组中块体单元的最大外接球半径为d3，d3的取值范围为：D/α2≤d3＜D/α3；4)以此类推，第n-1组中块体单元的最大外接球半径为dn-1，dn-1的取值范围为：D/αn-2≤dn-1＜D/αn-1；5)第n组中块体单元的最大外接球半径为dn，则dn的取值范围为：dn≥D/αn-1；其中，n＝int[log(D/d)/logα+0.416]+1,D为所有块体单元中最大块体单元外接球半径，d为所有块体单元中最小块体单元外接球半径，α＝2；步骤3.3，对第1组的块体单元与所有块体单元进行接触检测，具体步骤包括：1)取lmax＝2d1，将计算区域划分为以lmax为边长的正方体网格；2)根据公式：xK=1+int(xK-xminlmax+0.5)]]>yK=1+int(yK-yminlmax+0.5)]]>zK=1+int(zK-zminlmax+0.5)]]>将所有块体单元映射到网格上，其中，(xK,yK,zK)为第K个块体单元的形心坐标，K＝1,2,…,N，N为块体单元个数，xmin、ymin、zmin分别为计算区域在x、y、z方向的最小值；3)对每个块体单元与其周边所有相邻网格内块体单元进行接触检测，当循环到某个块体单元时，即以此块体单元所在网格为中心，判断此块体单元是否与该网格内部以及相邻网格内部的块体单元接触，其中，判断此块体单元是否与该网格内部以及相邻网格内部的块体单元接触的具体过程为：3-1)分别计算两个块体单元的形心之间的距离dist、两个块体单元各顶点与其对应的形心之间的最大距离，并计算两个块体单元各顶点与其对应的形心之间的最大距离的平均值lavmax；3-2)比较dist与lavmax，若dist≤2lavmax，则判断此两个块体单元接触，否则判断此两个块体单元不接触；步骤3.4，将第一组的块体单元从所有块体单元之中去除，对第2组块体单元与剩余第2组到第n组的块体单元进行接触检测，具体步骤包括：1)取lmax＝2d2，将计算区域划分为以lmax为边长的正方体网格；2)根据步骤3.3中2)至3)的方法，将第2组到第n组的块体单元映射到网格上，并对第2组块体单元与剩余第2组到第n组的块体单元进行接触检测；步骤3.5，以此类推，将已进行接触检测的块体单元组从所有块体单元中去除，并进行下一组块体单元与剩余所有块体单元组的接触检测，直至完成n组块体的接触检测。4.根据权利要求1所述的一种基于距离势函数三维任意凸多边形块体离散单元法，其特征在于，步骤四中，目标块体单元的法向接触力、切向接触力以及力矩具体计算步骤如下：步骤4.1，确定目标块体单元的最大内切球直径；步骤4.2，定义距离势函数，循环计算当前时间步每个与目标块体单元相接触的块体单元作用于目标块体的接触力和力矩；步骤4.3，计算所有与目标块体单元相接触的块体单元作用于目标块体的接触力和力矩。5.根据权利要求4所述的一种基于距离势函数三维任意凸多边形块体离散单元法，其特征在于，步骤4.2中，计算每个与目标块体相接触的块体单元作用于目标块体的接触力和力矩的具体步骤为：1)设块体单元βt与块体单元βc为两个相接触的块体单元，其中，块体单元βc为接触单元，βt为目标单元；2)定义距离势函数：其中，表示βt内部点p的势函数值，hp-1,hp-2,…,hp-M分别表示点p至βt各底面的距离，M为βt的底面数，h为βt的最大内切球半径；3)根据势函数法计算法向接触力，计算法向接触力，计算公式为：其中，fn是法向接触力；分别为βc、βt重叠区域内的点分别在βt与βc的距离势函数梯度；V为βt与βc的重叠区域；应用高斯公式，将法向接触力体...

【专利技术属性】
技术研发人员：赵兰浩，刘勋楠，毛佳，李同春，
申请(专利权)人：河海大学，
类型：发明
国别省市：江苏;32