基于离散单元法的提高颗粒离散接触检测效率的方法

摘要

本发明公开了属于颗粒离散元模型建立技术领域的一种基于离散单元法的提高颗粒离散接触检测效率的方法。该方法通过求解某时刻某颗粒最多接触的颗粒数，本方法中用到的数学模型都是一些简单的几何模型，通过对这些几何模型计算得到可接触颗粒数目与中心颗粒以及接触颗粒半径的关系，然后在程序中加入关系式来减少运算次数。可以使进行模拟计算的程序提前跳出循环，减少运算次数从而提高效率。使用该方法后计算量与现有模型相比有效减少，利用常用计算设备也可以进行大规模数值模拟计算。

说明书

技术领域

本发明属于颗粒离散元模型建立技术领域，特别涉及一种基于离散单元法的提高颗粒离散接触检测效率的方法。 

背景技术

离散单元法是由美国学者Cundall P.A.教授基于分子动力学原理首次提出的一种颗粒离散体物料分析方法，是一种显示求解的数值方法。自提出以来，在岩土工程和离散体工程这两大传统的应用领域中发挥了其他数值算法不可替代的作用。离散单元法是把介质看做由一系列离散的独立运动的单元组成的，根据离散物质本身所具有的离散特性建立数学模型，将需要分析的物体看做离散颗粒的集合。因此，离散单元法在分析具有离散体性质的物料时具有很大的优越性。但是离散单元法也存在着有别于其他数值方法的困难，比如，复杂的数据结构、分格检索、相邻颗粒的确定和颗粒接触产生或取消的判定等等，都将消耗大量微机内存和CPU时间，有时甚至会导致该用法的失败。这些都是实际工程应用中最基本的要求，也是离散单元法发展中需要解决的问题。 

在应用离散单元法进行数值模拟的过程中，把物料中的每个颗粒单独作为一个粒子单元建立数学模型，并给定粒子单元的尺寸和物理性质，如质量、刚度和阻尼等，各个粒子之间存在接触和分离这两种关系。当接触发生时，接触点处就会产生接触力和力矩，其大小可以根据接触力学模型求出。在模拟过程中，通过相邻颗粒之间的碰撞产生的接触力、力矩和不平衡力、力矩，计算每个颗粒在特定时刻的运动特征，通过每个颗粒的特征反映出整个系统的运动特性。离散单元法中的粒子单元具有一定的质量和形状，常用的散体颗粒几何形状是二维的圆形和三维的球形。 

模拟颗粒物质静态和运动行为时，需要根据颗粒间重叠量计算接触力，依此更新每个颗粒的速度和位置，进而确定性地演化整个颗粒体系，因此需要依据颗粒位置频繁计算颗粒是否发生重叠，计算强度巨大。目前常用的方法是划分网格， 网格法简洁明了，易于程序实现和并行化。然而接触检索关系的搜索仍然是颗粒离散元分析中最耗时间的部分。改进接触检测算法，提高检测效率是颗粒离散元模型的重要课题。 

发明内容

本发明的目的是提供一种基于离散单元法的提高颗粒离散接触检测效率的方法，其特征在于，包括以下几个部分： 

1)首先确立颗粒离散元模型；包括单一平面内的颗粒和立体颗粒，在循环流化床锅炉的气固两相流中的固体相颗粒的复杂流动过程还难以被人们充分了解，所以多采用数值模拟的方法来进行研究；将固体颗粒看作是离散的个体，来研究颗粒的运动，对每一个颗粒的受力情况和运动轨迹进行模拟；当考虑颗粒间碰撞时，将颗粒视为等直径的圆形或球体单元；当要计算颗粒i的运动轨迹时，先确定颗粒i的位置，再计算其他所有颗粒与颗粒j之间的距离L_ij；假设两颗粒i和j的其坐标分别为x_i，y_i和x_j，y_j表示，当两球心间的距离L_ij满足下式时就认为两球是接触的，发生了碰撞；其中D为颗粒直径； 

$>L_{\mathrm{ij}}=\sqrt{{(x_i-x_j)}^2+{(y_i-y_j)}^2}\le D,$>

2)所述单一平面内的颗粒，就是在二维情况下，颗粒离散在单一平面内，这里只涉及圆形颗粒；中心圆形颗粒与n个圆形颗粒外切；即简化成一个中心圆与两个外圆相切的情况： 

设中心圆半径为R，其他外切圆半径为r，连接三个圆心，中心圆圆心与两个外切圆圆心的连线夹角为θ，圆心角θ对应的弧长为l,其弧长由弧长公式得到，l＝θ(R+r)； 

根据勾股定理，得到下式： 

$>\cos \theta =\frac{{(R+r)}^2+{(R+r)}^2-{\left(2r\right)}^2}{2{(R+r)}^2}=1-\frac{{2r}^2}{{(R+r)}^2};$>

当有n个圆外切时，其总弧长为 

2π(R+r)＝nl 

$>n=\frac{2\pi }{\arccos [1-\frac{{2r}^2}{{(R+r)}^2}]}$>

实际情况下不会是正好n个圆与中心圆外切，所以在程序中令接触圆数目小于等于n即可； 

3)所述立体颗粒，即在三维情况下，颗粒离散在三个平面内，这里只涉及球形颗粒；连接球面上相切的三个球冠的圆心，构成一个非中心圆球球面上的球面三角形，设A为此三点与中心球球心所构成的正三棱锥的两个侧面所形成的二面角。则该球面正三角形的面积为S＝(R+r)²(3A-π)，其中三个二面角都相等；在这个正三棱锥中，过底面的一个顶点作对面侧棱的垂线，再连接另一个顶点；因为正三棱锥两个侧面三角形全等，与另一个顶点的连线必然也是垂线，所以这两个顶点间的底边垂直于该侧棱；这样就得到了正三棱锥二面角的平面角，即A；同时也是由于正三棱锥的两个侧面三角形全等，做出的这个三角形两个腰长相等，底边即为正三棱锥的底面三角形的一条边； 

下一步来求解角A，正三棱锥中侧面三角形的边分别为R+r、R+r和2r，则可求出θ，根据sin(θ/2)＝r/(R+r)，θ＝2arcsin[r/(R+r)]；所作的二面角A的平面角等腰三角形的腰长为(R+r)·sinθ，带入θ得： 

$>\left(\begin{array}{c}(R+r)\sin \theta =(R+r)\sin \left(2\arcsin \frac{r}{R+r}\right)=2(R+r)\sin \left(\arcsin \frac{r}{R+r}\right)\cos \left(\arcsin \frac{r}{R+r}\right)\\ =\frac{2r\sqrt{R^2+2\mathrm{Rr}}}{R+r}\end{array}\right)$>

底边就是该三棱锥与该棱异面的底边，2r。根据余弦定理来求A。 

$>\cos A=\frac{2{\left(\frac{2r\sqrt{R^2+2\mathrm{Rr}}}{R+r}\right)}^2-{\left(2r\right)}^2}{2{\left(\frac{2r\sqrt{R^2+2\mathrm{Rr}}}{R+r}\right)}^2}=\frac{R^2-r^2+2\mathrm{Rr}}{2R^2+4\mathrm{Rr}},A=\arccos \left(\frac{R^2-r^2+2\mathrm{Rr}}{2R^2+4\mathrm{Rr}}\right)$>

这个球面正三角形是由三个六分之一的球面圆(即球冠)与一个空缺组成的，意味着每个球面圆可以占据六份球面正三角形，而三个球面圆的六分之一在同一个球面正三角形里，而全球面的面积为4π(R+r)²，则全球面可以容纳的前述 球面正三角形数目为m＝4π(R+r)²/S，通过前面的分析可知，实际小球的数量大致是球面三角形的1/2(因边界处不能形成紧密的空缺)，故小球的数量为n＝0.5m(取不大于此结果的最大整数)； 

$>n=0.5m=\frac{0.5\times {4\mathrm{\pi R}}^2}{R^2\{3\arccos \left(\frac{R^2-r^2+2\mathrm{Rr}}{2R^2+4\mathrm{Rr}}\right)-\pi )\}}=\frac{2\pi }{3\arccos \left(\frac{R^2-r^2+2\mathrm{Rr}}{2R^2+4\mathrm{Rr}}\right)-\pi }$>

经过上述计算，得到了三维情况下颗粒最多可接触的颗粒，这样，有了二维和三维情况最多可接触的颗粒数目，在对循环流化床中固相颗粒进行数值模拟计算时，就可以将得到的结果写进程序来优化计算；而且应用范围不仅仅局限在循环流化床中，在其他需要判断颗粒碰撞的工程中都可以使用。 

本发明的有益效果是通过求解某时刻某颗粒最多接触的颗粒数，通过简单的几何模型来判断某颗粒(圆或球)周围最多可接触的颗粒，可以使进行模拟计算的程序提前跳出循环，减少运算次数从而提高效率，使大规模的数值模拟节省计算时间。 

附图说明

图1是二维情况下的圆形颗粒，选取任意两个外切圆来与中心圆组成该模型。其中，R为中心圆的半径，r为外切圆的半径，θ为圆心角，l为θ对应的弧长。 

图2是三维情况下的球形颗粒，选取任意三个外切球来与中心球组成该模型。其中，R为中心球的半径，r为外切球的半径。 

图3是图二中连接四个球心后形成的正三棱锥，R+r为正三棱锥的棱长，2r为正三棱锥底面正三角形的边长，θ为正三棱锥侧面三角形的顶角，A为两个侧面形成的两面角。 

图4是对循环流化床气固两相流的模拟图，模拟对象即为里面的颗粒。 

具体实施方式

下面结合附图对本发明做进一步说明。 

图1-3所示为一种基于离散单元法的提高颗粒离散接触检测效率的方法，通过求解某时刻某颗粒最多接触的颗粒数，过程就是依靠简单的几何模型来判断某 颗粒(圆或球)周围最多可接触的颗粒，模拟颗粒静止和运动行为时，需要根据颗粒间重叠量计算接触力进而求出每个颗粒的速度和位置，频繁计算颗粒是否发生重叠，其过程为： 

(1)二维情况(圆形颗粒) 

对于二维颗粒，检索某颗粒与其他颗粒接触时，传统的模拟计算是对每两个颗粒都要进行一次判断，所以通过计算某颗粒最多能够接触的颗粒数目，来提前结束程序循环。 

如图1所示，三个圆互相相切。左侧的圆为中心圆，半径为R。右侧的圆是外切圆，半径为r。θ角是连接三个圆心，形成的左侧圆的圆心角。L是θ角对应的弧长。 

第一步利用勾股定理来计算圆心角θ 

$>\cos \theta =\frac{{(R+r)}^2+{(R+r)}^2-{\left(2r\right)}^2}{2{(R+r)}^2}=1-\frac{{2r}^2}{{(R+r)}^2}$> ① 

第二步利用弧长公式来计算圆心角θ对应的弧长 

l＝θ(R+r) ② 

第三步假设n为外切圆的数目(先假设正好有n个圆与中心园外切)，则 

2π(R+r)＝nl ③ 

连立上面三个公式，求出n 

$>n=\frac{2\pi }{\arccos [1-\frac{{2r}^2}{{(R+r)}^2}]}$>

求出的n即为某圆形颗粒最多可以接触到的颗粒。考虑到一般情况下不会是正好n个圆与中心圆外切，所以可以在程序中取不大于n的最大整数来减少运算。 

(2)三维情况(球形颗粒) 

对于三维颗粒，检索某颗粒与其他颗粒接触时，传统的模拟计算是对每两个颗粒都要进行一次判断，所以通过计算某颗粒最多能够接触的颗粒数目，来提前结束程序循环。 

如图2、图3所示四个球互相相切。下面的球为中心球，半径为R。上面的 球为外切球，半径为r。连接四个球心。形成一个正三棱锥，见图三。三条棱长为R+r，底面等边三角形边长为2r。θ角是正三棱锥的侧面的一个顶角。A角是两个侧面形成的两面角。 

球冠：球面被平面所截得的一部分叫做球冠.截得的圆叫做球冠的底，垂直于截面的直径被截得的一段叫做球冠的高。球冠也可以看作一段圆弧绕经过它的一个端点的直径旋转所成的曲面。 

球面三角形：把球面上的三个点用三个大圆弧联结起来，所围成的图形叫做球面三角形。面积公式：S＝R²(A+B+C-π)，其中R为球的半径，A、B、C分别是球面三角形的三个两面角。 

二面角：平面内的一条直线把平面分为两部分，其中的每一部分都叫做半平面，从一条直线出发的两个半平面所组成的图形，叫做二面角(这条直线叫做二面角的棱，每个半平面叫做二面角的面)。 

1)连接球面上相切的三个球冠的圆心，构成一个球面三角形(此时该球面已经不是中心球上的球面了，是以R+r为半径构成的大球)。A为此三点与中心球球心所构成的正三棱锥的两个侧面所形成的二面角。则该球面正三角形的面积为S＝(R+r)²(3A-π)(三个二面角都相等)。在这个正三棱锥中，过底面的一个顶点作对面侧棱的垂线，再连接另一个顶点。因为正三棱锥两个侧面三角形全等，与另一个顶点的连线必然也是垂线，所以这两个顶点间的底边垂直于该侧棱；这样就得到了正三棱锥二面角的平面角，即A；同时也是由于正三棱锥的两个侧面三角形全等，做出的这个三角形两个腰长相等，底边即为正三棱锥的底面三角形的一条边。 

2)求解角A。正三棱锥中三个的侧面三角形的两个腰长为R+r，底边为2r，则根据正弦公式可以求出θ。 

sin(θ/2)＝r/(R+r) 

θ＝2arcsin[r/(R+r)] 

因为二面角A的平面角等腰三角形的腰长为(R+r)·sinθ，带入θ得： 

$>\left(\begin{array}{c}(R+r)\sin \theta =(R+r)\sin \left(2\arcsin \frac{r}{R+r}\right)=2(R+r)\sin \left(\arcsin \frac{r}{R+r}\right)\cos \left(\arcsin \frac{r}{R+r}\right)\\ =\frac{2r\sqrt{R^2+2\mathrm{Rr}}}{R+r}\end{array}\right)$>

同时二面角A的平面角等腰三角形的底边为2r。根据余弦定理来求A。 

$>\cos A=\frac{2{\left(\frac{2r\sqrt{R^2+2\mathrm{Rr}}}{R+r}\right)}^2-{\left(2r\right)}^2}{2{\left(\frac{2r\sqrt{R^2+2\mathrm{Rr}}}{R+r}\right)}^2}=\frac{R^2-r^2+2\mathrm{Rr}}{2R^2+4\mathrm{Rr}},A=\arccos \left(\frac{R^2-r^2+2\mathrm{Rr}}{2R^2+4\mathrm{Rr}}\right)$>

3)求解数目。做出的球面正三角形是由三个六分之一的球面圆(即球冠)与一个空缺组成的。意味着每个球面圆可以占据六份球面正三角形，而三个球面圆的六分之一在同一个球面正三角形里。而全球面的面积为4π(R+r)²，则全球面可以容纳的前述球面正三角形数目为m＝4π(R+r)²/S。 

$>n=0.5m=\frac{0.5\times {4\mathrm{\pi R}}^2}{R^2\{3\arccos \left(\frac{R^2-r^2+2\mathrm{Rr}}{2R^2+4\mathrm{Rr}}\right)-\pi )\}}=\frac{2\pi }{3\arccos \left(\frac{R^2-r^2+2\mathrm{Rr}}{2R^2+4\mathrm{Rr}}\right)-\pi }$>

通过前面的分析可知，实际小球的数量大致是球面三角形的1/2(因边界处不能形成紧密的空缺)，故小球的数量为n＝0.5m，所以可以在程序中取接触球数目为不大于n的最大整数来减少运算。 

在对循环流化床中的固相颗粒进行数值模拟时，使用的计算机处理器为Inter(R)Core(TM)i5-2400,主频为3.10GHz；RAM为4GB。使用的软件为Fortran。模拟对象为150mm×4mm×900mm矩形截面的准三维流化床(床深度为颗粒直径)，网格划分为17×92，颗粒采用2400个直径均为4mm的小球。没考虑可接触的颗粒数目，单纯进行循环计算，对每两个颗粒都进行判断时，得到颗粒的速度分布用时为25小时40分钟；经过该发明优化程序后，得到同样的结果只耗时20小时18分钟。效率相比之前提高了20.91％。可见，经过改发明的优化后，效率提高是比较明显的。这样，没有通过提高计算机硬件以及并行计算，单独通过修改程序，就达到了提高计算效率的目的。