一种基于离散元法的弯管中颗粒流动问题的数值模拟

刘 红

(安徽理工大学理学院，安徽 淮南 232001)

一种基于离散元法的弯管中颗粒流动问题的数值模拟

刘 红

(安徽理工大学理学院，安徽 淮南 232001)

本文建立了直角弯管中颗粒流动问题的数学模型，并对核心模块的算法进行详细描述，在此基础上设计了基于离散元法求解弯管中颗粒流动问题的算法与程序。数值实验结果表明：实验结果与实际物理现象相符，验证了所设计的求解算法和所研制的程序的正确性。

离散元法；直角弯管；颗粒

1 引言

管道运输具有一些其他运输机械无法比拟的优势,它在工农业和物流运输业中应用非常广泛[1],目前在高速运输的数值模拟方面和管壁磨损方面的研究是较多的, 对于颗粒在运动过程中对弯管的磨损机理研究也逐渐被人们重视[2]。弯管中颗粒的流动增加了问题处理的复杂度,因此多数选取球形颗粒(二维选取圆盘形颗粒)来模拟。

近年来,计算机软件和硬件的快速发展[3],为长时间大规模的数值模拟打下了坚实的基础。离散单元法在研究气-固两相流中固相方面有独特的优势, 离散元这一概念最早是1971[4]年Cundall为了分析和处理岩土工程方面的问题而提出来的,对于不连续介质，离散元法是行之有效的数值模拟方法[5]。在实际应用方面,离散元技术在岩土、化工、农业、风沙流动及散粒材料的运输等领域有广泛应用。1986年王泳嘉[6]首次在会议上向国内岩石力学与工程界介绍了离散元法。李世海[7]在2012年颗粒材料计算力学会议上发表了CDEM的最新研究进展。Zhang[8-9]等对流体输送过程用离散元方法进行了数值模拟,分析了颗粒之间的相互作用和颗粒对管道的冲蚀磨损作用,并侧重预测输送管道最容易损坏的位置,对实际的流体输送工程问题有很大的参考价值。Yue等[10]基于CPU-GPU异构架构将DEM的串行代码做了并行化处理，提升了其的计算效率。毛靖儒等[11]采用离散元法对弯管内稀疏气固两相流进行了模拟，同时他也研究了固体颗粒对壁面磨损效应,指出了壁面处粒子的分布、碰撞角度及浓度分布是固粒对壁面磨损的主要因素。

本文首先对Hertz模型和M-D模型进行详细的总结并给出了平动、转动平衡方程的离散格式，然后给出模型的基本假设，并建立直角弯管中颗粒流动问题的模型，给出了基于离散元法的核心模块的算法,最后用直角弯管中颗粒流动的数值实验对我们的程序设计进行了正确性验证，得出结论。

2 离散单元法

考虑到运动颗粒之间的碰撞并使用相互作用方程，颗粒受到系统中其它颗粒对其的接触作用,遵循质心运动定理且满足牛顿第二运动定律，则有

(1)

(2)

(3)

其中E1(E2)和,μ1(μ2)分别表示颗粒的弹性模量和泊松比。

粒间的法向接触力增量ΔF与两颗粒之间的重叠量Δα关系如下:

(4)

在实际的接触问题中,颗粒间存在着相对滑动和摩擦作用,且这些作用是不可以忽略的,因此Hertz模型不能够准确地描述接触力,由于相互接触产生的切向接触力的作用,使得两颗粒在切向力方向上发生相对滑移,滑移从接触表面开始之后蔓延到颗粒的内部,Mindlin-Deresiewic模型考虑了切向力与切向位移的联系,可以更加真实地描述这一过程。切向力增量为

(5)

其中，δ，Δδ和T分别表示两个颗粒接触表面切向位移，切向位移增量和切向力。通常采用向前差分(中心差分)格式对方程左端的导数项进行离散，用显示格式对平动和转动平衡方程的进行求解 ,因此需要选取合适的时间步长使得两颗粒间的作用力在Δt的时间内变化非常小，这样计算的稳定性和准确性才能得到保证。通常我们由Rayleigh波的波速来确定时间步长。假设颗粒的材质为弹性固体, ρ,G，ν，R分别表示颗粒的密度，剪切模量、泊松比、不同尺寸的颗粒组成的系统,时间步长由下式计算得出

(6)

式(6)是在颗粒间相对速度较小时计算得出的,因此这种确定时间步长的方法比较适用于颗粒运动不太剧烈的系统。依据实际情况选取合适的时间步长来确保稳定性,可以在时间步长前乘一个常数C (0

3 直角弯管中颗粒流动模型的建立与算法描述

3.1 物理模型的建立和基本参数

二维平面内的直角弯管,实际上也就是同心圆环的四分之一圆周,为此我们将圆环进行剖分,采用一小段线段代替圆周上所对应的一段圆弧。根据二维平面圆的参数方程(7)即可得到每一面墙的两个端点坐标。

(7)

用球形颗粒代表真实颗粒,用墙表示管道壁,建立离散元模型模拟了直角弯管中颗粒在重力作用下的流动,表1为本次数值模拟的参数，图1表示初始状态.

表1 颗粒和墙的物理参数

图1 初始状态图

对弯管颗粒流动的离散元模拟,首先要制备基本稳定的颗粒体系。在颗粒生成的团块区域随机的生成1000个颗粒,并让颗粒在重力作用下自然沉降,达到稳定状态,如图1所示，并将该时视为零时刻。下面主要介绍求解过程中最重要的颗粒运动模块，接触判断模块和接触力求解模块。

运动模块的主要功能为计算颗粒在每一时刻的速度与位移。从tn时刻到tn+1时刻对颗粒i的运动情况可由算法1描述。

算法1：

颗粒间接触力求解模块利用前面介绍的接触模型计算颗粒的法向接触力和切向接触力。算法2为求解颗粒i与颗粒j之间接触力的算法。

高同型半胱氨酸（Hcy）是半胱氨酸和蛋氨酸代谢的一种重要中间产物，是一种含硫非必需氨基酸，能够损伤血管内皮细胞，是动脉粥样硬化、血栓、糖尿病肾病的独立危险因素[7-8]。Hcy由肾脏合成与代谢，当肾小球滤过功能、肾小管代谢功能损害时，血液Hcy水平升高，可以与高血糖、高糖化血红蛋白协同损伤肾脏血管内皮细胞、肾小球基底膜细胞，使肾小球滤过膜孔径增大，造成肾滤过功能损害，加重蛋白尿程度。

算法2：

Step 1求颗粒之间的重叠量α：

2) 计算两球形颗粒间的重叠量α=ri+rj-d,式中rl为颗粒l(l=i，j)的半径。

Step 2 计算法向接触力fn+1、法向力Fn+1、接触圆半径an+1以及法向接触力增量ΔFn+1。

Step 4求合力以及合力矩：

颗粒之间的法向力是有Hertz模型求得的,其算法描述如下。

算法3：

Step 3 生成法向力大小：Fn+1=fn+1+pn+1。

颗粒间切向力大小是由M-D模型计算出来的,其算法描述如下。

算法 4：

Step 3 切向接触力大小：tn+1=tn+Δtn+1。

4 直角弯管中颗粒流动数值模拟

打开挡板,颗粒在重力作用下,沿着弯管向下流动,并在最终出口向外流出, 图2给出流动过程中五个不同时刻状态图,其中t1=0.23842s，t2=0.47684s，t3=0.71526s，t4=0.95368s，t5=1.1921s。

图2 五个不同时刻的颗粒流动状态

如图2所示,在初始状态,颗粒由于重力的作用堆积在挡板上面。打开挡板之后,颗粒在重力作用下首先向下落,由于在下落过程中颗粒之间会相互碰撞，进而产生接触力使得颗粒的运动速度发生变化,在这一过程中有靠近管壁面的部分颗粒也会与管壁发生接触,壁面附近有少量颗粒集中,这就导致了不同颗粒有不同的运动状态,所以在竖直管下落的阶段,从整体上看颗粒由集密变得稀疏。

当颗粒流前端接触到弯管,颗粒与颗粒、颗粒与弯管壁之间发生碰撞并出现聚团现象。从图2的第3个和第4个图可以看出,颗粒流进入弯管后,由于向心力作用大部分颗粒沿着弧形弯管外壁运动；颗粒从竖直管进入弯管的过程中,颗粒在与壁面发生碰撞之前其运动的方向改变较小,颗粒与壁面接触后并在向心力的作用下在弯管外壁附近移动,当然也有少部分颗粒受到弯管外壁的弹力作用往弯管内壁方向运动,这一过程中有的颗粒也会与弯管内壁发生接触碰撞；图2上很少颗粒在弯管内壁附近运动。上述结果说明弯管外管壁承是颗粒流的冲击和磨损作用的主要作用区域。

颗粒离开弯管后进入水平管道,由于重力的原因颗粒都集中在下管壁流动。可以看出颗粒间的碰撞基本恢复均匀。最后颗粒离开管道做抛物线运动,落在水平地面上，水平速度的不同使得它们的抛物线形状也有所不同。

5 结 论

本文首先对弯管结构进行了数学建模,并用离散元法模拟了直角弯管中颗粒在重力作用下的流动现象，最后的实验结果表明设计的程序能够准确地模拟弯管中颗粒流动这一模型。

[1] Stoss HA．Pneumatic Conveying[M]．New York：Wiley Interscience，1983．

[2] 卢洲，刘雪东，潘兵．基于CFD-DEM方法的柱状颗粒在弯道中输送过程的数值模拟[J]．中国粉体技术，2011，17(5)：1-2．

[3] 石钟慈．第三种科学方法—计算机时代的科学计算[M]．清华大学出版社，2000．

[4] Cundall P A．A computer model for simulating Progressive large scale movements in blocky R- rock systems[C]// Muller Led．Proc Symp Int Soc Rock Mechanics．Rotterdam：Balkama A A，1971 (1)：1-8．

[5] 王泳嘉，邢纪波．离散单元法同拉格朗日元法及其在岩土力学中的应用[J]．岩土力学，1979，16(2)：1-6．

[6] 王泳嘉．离散单元法-一种适用于节理岩石力学分析的数值方法[A]//第一届全国岩石力学数值计算机模型实验讨论会会议论文集[C]，1986：32-37．

[7] 李世海，冯春，刘晓宇，等．用于描述地质体块体颗粒的连续及非连续计算模型-CDEM最新进展[A]．颗粒材料计算力学研究进展会议录[C]，2012．

[8] Hao Zhang，Yuanqiang Tan，Dongmin Yang，et al． Numerical investigation of the location of maximum erosive wear damage in elbow：Effect of slurry velocity，bend orientation and angle of elbow[J]．Technology．2012，217(3)：467- 476．

[9] Yuanqiang Tan，Hao Zhang．Numerical simulation of concrete pumping process and investigation of wear mechanism of the piping wall[J]．Tribology international．2012，46(1)：137-144．

[10] Xiaoqiang Yue，Hao Zhang，Congshu Luo，et al．Parallelization of a DEM code Based on CPU-GPD Heterogeneous Architecture[J]．Paraller Computational Fluid Dynamics，2014，405：149-159．

[11] 柳成文，毛靖儒，俞茂铮．90°弯管内稀疏气固两相流及固粒对壁面磨损量的数值研究[J]．西安交通大学学报，1999，33(9)：53-57．

[12] 孙其诚，王光谦．颗粒物质力学导论[M]．北京：科学出版社，2009．

[责任编辑：崔海瑛]

(英文摘要略)

DOI 10.13356/j.cnki.jdnu.2095-0063.2016.06.020

Numerical simulation of particle flow in bend pipe based on discrete element method

LIU Hong

(School of Science in Anhui University of Science and Technology, Huainan 232001, China)

刘红(1992-)，女，河南信阳人，硕士研究生，主要从事Petri网理论与应用及并行计算方面的研究。

TQ022.3

A

2095-0063(2016)06-0089-05

2016-06-28

DOI 10.13356/j.cnki.jdnu.2095-0063.2016.06.019