上海理工大学学报  2023, Vol. 45 Issue (6): 543-551, 601   PDF    
非规则颗粒与滚筒球磨机相互作用的离散元–有限元耦合分析
李天泽, 王嗣强, 杨冬宝, 季顺迎     
大连理工大学 工业装备结构分析优化与CAE软件全国重点实验室,大连 116024
摘要: 非规则颗粒形态显著影响球磨机中颗粒材料的运动行为,同时对球磨机结构产生不同的力学响应。采用超二次曲面方程描述非规则颗粒材料,并发展了离散元(discrete element method, DEM)–有限元(finite element method, FEM)耦合算法以分析非规则颗粒材料与滚筒球磨机间的相互作用。在该算法中,滚筒球磨机被划分为一系列四边形单元,并且超二次曲面颗粒与球磨机间的接触判断可转化为超二次曲面颗粒与四边形单元间的接触判断。将计算得到的碰撞力向节点插值,该节点力作为载荷条件在有限元中进行隐式动力学求解,并通过求解颗粒运动方程和有限元动力方程实现颗粒运动更新和结构力学响应分析。进一步分析了滚筒球磨机中球体、圆柱体和立方体颗粒材料的运动行为以及滚筒球磨机的变形及应力分布规律。研究结果表明:颗粒形状显著影响球磨机研磨效果和球磨机结构位移;立方体颗粒具有最好的研磨效果,而圆柱体颗粒比球体颗粒具有更好的研磨效果;在结构位移上,立方体颗粒有最大的峰值位移,其次是圆柱体和球体颗粒;相比于规则的球体颗粒,超二次曲面颗粒有不规则的表面形状,这会导致更好的研磨效果和更大的结构位移。
关键词: 离散元方法     有限元方法     耦合算法     超二次曲面单元     非规则颗粒    
A DEM-FEM coupled analysis for the interaction between the irregular granular materials and tumbling mill
LI Tianze, WANG Siqiang, YANG Dongbao, JI Shunying     
State Key Laboratory of Structural Analysis, Optimization and CAE Software for Industrial Equipment, Dalian University of Technology, Dalian 116024, China
Abstract: The irregular granular morphology significantly affects the kinematic behavior of the granular materials in a tumbling mill, and simultaneously causes different mechanical responses of the tumbling mill. The superquadric functions were used to describe the irregular granular materials, and a discrete element method (DEM)-finite element method (FEM) coupled algorithm was developed to analyze the interaction between the irregular granular materials and tumbling mill. In this algorithm, the tumbling mill was divided into a series of quadrilateral elements, and the contact detection between the superquadric particles and the tumbling mill could be simplified as the contact detection between the superquadric particles and the quadrilateral elements. The calculated collision forces were interpolated in the nodes, which were used as the loading conditions in the finite elements for the solution of the implicit dynamics. The update of particle motion and the structural mechanical response were realized by solving the particle motion equation and the finite element dynamic equation. Furthermore, the motion behavior of spherical, cylindrical and cubic granular materials in a tumbling mill and the deformation and stress distribution of the tumbling mill were analyzed. The numerical results indicated that the particle shape significantly affected the grinding effect and structural displacement of the ball mill. The cube particles had the best grinding effect, while cylindrical particles had better grinding effect than spherical particles. In terms of structural displacement, cubic particles had the largest peak displacement, followed by the cylindrical and spherical particles. Compared to regular spherical particles, superquadric particles had irregular surface shapes, which lead to a better grinding effect and greater structural displacement.
Key words: discrete element method     finite element method     coupled algorithm     superquadric element     irregular granular materials    

滚筒球磨机是采矿工程中以颗粒材料为研磨介质的通用研磨设备,因此,内部颗粒的几何形态显著影响其研磨性能[1]。一般采用规则的球体颗粒材料作为研磨介质,但根据不同工程应用需求,非规则形态的颗粒也被广泛使用,如鹅卵石形状[2]、矿石形状[3]、圆柱体[4]和立方体[5]等。即便是采用球体颗粒,在滚筒球磨机的旋转过程中,球体颗粒材料也会经历从最初规则形状到非规则颗粒混合物的形态转变。这主要是由于碰撞和多种磨损机制导致规则球体颗粒瞬间磨损成非规则形态颗粒[6-7]。同时,在滚筒球磨机中球体和非规则颗粒材料具有不同的运动特性,并且从球体颗粒材料得到的结论难以直接应用于非规则颗粒材料[8-11]。因此,探究滚筒球磨机内非规则颗粒材料的运动规律及结构的应力分布规律,对于球磨机的设计和优化至关重要。

自1979年Cundall等[12]首次提出离散单元法(discrete element method, DEM)以来,经过40多年的发展,它已经成为了模拟颗粒材料力学特性的重要工具[13-15]。在滚筒球磨机的离散元模拟中,大量球体颗粒的材料参数及运动状态被广泛研究,这对于球磨机的参数选择和结构优化具有重要的应用价值。为了更准确地模拟颗粒材料的运动状态,非规则颗粒材料的离散元方法迅速发展。模拟非规则颗粒主要采用两种方法:单个非规则颗粒模型和组合颗粒模型。单个非规则颗粒模型使用一个连续包络函数或几何拓扑方法来构造非规则颗粒,如椭球体[16]、超二次曲面[17]和扩展多面体[18];组合颗粒模型则采用多个凸形颗粒以任意空间位置来形成具有复杂形态的非规则颗粒。在这些方法中,基于连续函数包络的超二次曲面单元可以合理地模拟自然界中约80%的颗粒形状[19]。通过调整函数中的半轴长和形状参数,可以改变模拟的颗粒尺寸和表面尖锐度[20-22]。这种方法通过联立求解方程组进而计算颗粒间和颗粒与结构间的接触点,并通过非线性接触力模型计算碰撞力。然而,在研磨过程中球磨机结构的变形和应力等力学响应对球磨机的安全运行至关重要。工程结构分析中,有限单元法(finite element method, FEM)是一种广泛应用的数值方法,其基于宏观连续介质力学模型并通过建立材料的本构关系来描述物理特性,随后求解控制方程以获得结构的变形和应力分布等[23]。因此,为了分析球磨机的变形和应力分布特性,有必要采用有限单元法对滚筒球磨机进行建模和数值分析。

值得注意的是,采用单一离散元方法或有限元方法难以准确模拟非规则颗粒材料在滚筒球磨机中的运动、碰撞以及与球磨机的相互作用过程。因此,本文采用超二次曲面方程构造非规则颗粒形态,建立滚筒球磨机的有限元结构模型,同时发展了DEM-FEM耦合算法,以分析不同形状颗粒在球磨机中的运动行为,以及滚筒球磨机的内部结构变形及应力分布规律。

1 超二次曲面DEM-FEM耦合算法

采用超二次曲面方程构造球体和非规则形态颗粒,同时发展超二次曲面颗粒与四边形面间的搜索算法以计算接触点和重叠量,并通过非线性接触力模型计算颗粒间及颗粒与结构间的碰撞力。最后,通过求解有限元动力方程及颗粒运动方程进而实现结构位移–应力计算和颗粒运动更新。

1.1 基于超二次曲面方程的非规则颗粒构造

超二次曲面方程是一种椭球方程的扩展形式,可通过改变其5个参数进而得到颗粒的不同尺寸和形态。以颗粒中心为原点,其函数形式可表示为[24-25]

$ {f\left(x,y,z\right)=\left({\left|\frac{x}{a}\right|}^{{n}_{2}}+{\left|\frac{y}{b}\right|}^{{n}_{2}}\right)}^{{n}_{1}/{n}_{2}}+{\left|\frac{z}{c}\right|}^{{n}_{1}}-1=0 $ (1)

式中:$ a,b $$ c $分别表示颗粒沿3个主轴方向的半轴长;$ {n}_{1} $$ {n}_{2} $为形状参数。

图1为采用不同半轴长和形状参数得到的不同形态颗粒。在形状参数的选择上:当$ {n}_{1}= {n}_{2}= 2 $时,可以生成球体或椭球体;当$ {n}_{1}\gg 2 $$ {n}_{2}=2 $时,可以生成类圆柱体;当$ {n}_{1}={n}_{2}\gg 2 $时,可以生成类立方体。对于超二次曲面颗粒间接触点的计算方法可参考文献[17],并且该牛顿迭代方法具有较好的数值精度。


图 1 由超二次曲面方程构造的颗粒形状 Fig. 1 Particle shapes constructed by the superquadric equations
1.2 超二次曲面颗粒与结构间的接触检测算法

在离散元–有限元耦合模拟中,不同形态颗粒与结构间的参数传递是至关重要的。首先,通过颗粒与结构间的接触检测算法计算接触点和重叠量,同时采用非线性接触力模型计算碰撞力和力矩;然后,采用形函数将碰撞力插值至有限单元的节点上;最后,通过求解有限单元的动力方程及颗粒运动方程进而更新颗粒运动信息和结构信息。超二次曲面颗粒与复杂结构间的接触检测可转化为超二次曲面颗粒与每个实体单元的四边形面间的接触判断,同时这种简化后的接触判断存在3种接触模式:面接触、边接触和点接触,如图2所示。图中:${{\boldsymbol{a}}},{{\boldsymbol{b}}},{{\boldsymbol{c}}},{{\boldsymbol{d}}}$分别表示点B坐标减点A坐标、点 C坐标减点 B 坐标、点 D 坐标减点 C 坐标、点 A坐标减点 D 坐标得到的坐标向量;${{{\boldsymbol{n}}}}_{\mathrm{w}}$为平面的单位法向量,${{{\boldsymbol{n}}}}_{\mathrm{w}}= {{\boldsymbol{a}}}\times \dfrac{\boldsymbol{b}}{\left|{\boldsymbol{a}}\times {\boldsymbol{b}}\right|}$;点$ A $指向超二次曲面中心点$ P $的向量可表示为${{\boldsymbol{e}}}$,颗粒中心到平面$ \left|P{Q}_{P}\right| $的距离可表示为$\left|{{\boldsymbol{e}}}{\text{·}} {{{\boldsymbol{n}}}}_{\mathrm{w}}\right|$


图 2 超二次曲面单元与四边形面的接触判断 Fig. 2 Contact detection between a superquadric particle and a quadrilateral element

超二次曲面颗粒与四边形面之间的边接触和点接触算法与文献[17]中所介绍的接触算法一致。本文将重点介绍超二次曲面颗粒与四边形面之间的接触算法。首先,对超二次曲面颗粒的包围球与四边形单元所在平面进行粗判断。如果包围球发生接触,则需要计算超二次曲面颗粒与平面间的精确接触点,并判断该接触点是否位于四边形面内。一个平面方程可表示为

$ {\boldsymbol{n}}_{\mathrm{w}}{\text{·}} ({\boldsymbol{x}}-{\boldsymbol{x}}_{A})=0 $ (2)

式中,${\boldsymbol{x}}_{A}$为平面上点$ A $的坐标。

超二次曲面颗粒的包围球半径设定为$ {r}_{\mathrm{p}} $,则采用包围球方法的粗判断可表示为

$ \left\{\begin{array}{l}{r}_{\mathrm{p}}-\left|{{\boldsymbol{e}}}{\text{·}} {{{\boldsymbol{n}}}}_{\mathrm{w}}\right|\leqslant 0,接触\\ {r}_{\mathrm{p}}-\left|{{\boldsymbol{e}}}{\text{·}} {{{\boldsymbol{n}}}}_{\mathrm{w}}\right| > 0,不接触\end{array}\right. $ (3)

如果超二次曲面颗粒的包围球与平面发生接触,则需要确定超二次曲面颗粒与平面的接触点位置。假定颗粒表面上一点${{\boldsymbol{x}}}{\left(x,y,z\right)}^{\mathrm{T}}$与平面间的距离取最大值,其优化方程表示为[26]

$ \left\{\begin{array}{l}{\mathrm{max}} ({{{\boldsymbol{n}}}}_{\mathrm{w}}{\text{·}} {\boldsymbol{x}})\\ F\left({{\boldsymbol{x}}}\right)=0\end{array}\right. $ (4)

式中:$F\left({{\boldsymbol{x}}}\right)=f\left[{{{\boldsymbol{R}}}}^{\mathrm{T}}\cdot \left({{\boldsymbol{x}}}-{{{\boldsymbol{x}}}}_{\mathrm{c}}\right)\right]$为全局坐标系下超二次曲面函数方程;${{\boldsymbol{x}}}_{\mathrm{c}}$为超二次曲面单元的质心位置;${{\boldsymbol{R}}}$为基于四元数的转换矩阵。取中间变量$ \tilde{x} $,通过代入$ x=a\tilde{x} $$ y={\alpha }_{1}b\tilde{x} $$ z={\beta }_{1}c\tilde{x} $$ {\gamma }_{1}={\left(1+{\alpha }_{1}^{{n}_{2}}\right)}^{{n}_{1}/{n}_{2}-1} $,并考虑平面法向量${\boldsymbol{n}}_{\mathrm{w}}\left({n}_{x},{n}_{y},{n}_{z}\right)$的分量具有正负符号,可将超二次曲面颗粒与平面的解析解表示为

$ \left\{\begin{array}{l}{\alpha }_{1}={\left(\left|b{n}_{y}\right|/\left|a{n}_{x}\right|\right)}^{1/\left({n}_{2}-1\right)}\\ {\gamma }_{1}={\left(1+{\alpha }_{1}^{{n}_{2}}\right)}^{{n}_{1}/{n}_{2}-1}\\ {\beta }_{1}={\left({\gamma }_{1}\left|{n}_{z}c\right|/\left|{n}_{x}a\right|\right)}^{1/\left({n}_{1}-1\right)}\\ x=a\Big/{\left[{\left(1+{\alpha }_{1}^{{n}_{2}}\right)}^{{n}_{1}/{n}_{2}}+{\beta }_{1}^{{n}_{1}}\right]}^{1/{n}_{1}}{\text{·}} {\rm{sign}}\left({n}_{x}\right)\\ y={\alpha }_{1}b\left|x\right|/a{\text{·}} {\rm{sign}}\left({n}_{y}\right)\\ z={\beta }_{1}c\left|x\right|/a{\text{·}} {\rm{sign}}\left({n}_{z}\right)\end{array}\right. $ (5)

如果满足$\left[\left({\boldsymbol{x}}_{A}-{{\boldsymbol{x}}}\right){\text{·}} {\boldsymbol{n}}_{\mathrm{w}}\right]\leqslant 0$,则颗粒与四边形单元所在的平面发生接触,且超二次曲面颗粒与平面的重叠量可表示为

$ {{{\boldsymbol{\delta}} }}_{\mathrm{n}}=\left[\left({{{\boldsymbol{x}}}}_{A}-{{\boldsymbol{x}}}\right){\text{·}} {{{\boldsymbol{n}}}}_{\mathrm{w}}\right]{\text{·}} {{{\boldsymbol{n}}}}_{\mathrm{w}} $ (6)

此外,为了进一步判断接触点是否在四边形面内,这里取4个未知参数$ {\phi }_{1},{\phi }_{2},{\phi }_{3} $$ {\phi }_{4} $,可分别表示为

$ \left\{\begin{array}{l} {\phi}_{1}=\left[\left({\boldsymbol{x}}+{{\boldsymbol{\delta }}}_{\mathrm{n}}-{{\boldsymbol{x}}}_{A}\right)\times {{\boldsymbol{a}}}\right]\cdot {{\boldsymbol{n}}}_{\mathrm{w}}\\ {\phi }_{2}=\left[\left({\boldsymbol{x}}+{{\boldsymbol{\delta }}}_{\mathrm{n}}-{{\boldsymbol{x}}}_{B}\right)\times {{\boldsymbol{b}}}\right]\cdot {{{\boldsymbol{n}}}}_{\mathrm{w}}\\ {\phi }_{3}=\left[\left({\boldsymbol{x}}+{{\boldsymbol{\delta }}}_{\mathrm{n}}-{{\boldsymbol{x}}}_{C}\right)\times {{\boldsymbol{c}}}\right]\cdot {{{\boldsymbol{n}}}}_{\mathrm{w}}\\ {\phi }_{4}=\left[\left({\boldsymbol{x}}+{{\boldsymbol{\delta }}}_{\mathrm{n}}-{{\boldsymbol{x}}}_{D}\right)\times {{\boldsymbol{d}}}\right]\cdot {{{\boldsymbol{n}}}}_{\mathrm{w}} \end{array}\right.$ (7)

式中:$ {{\boldsymbol{x}}}_{B},{{\boldsymbol{x}}}_{C} $为和$ {{\boldsymbol{x}}}_{D} $分别为平面上点$ B $$ C $$ D $的坐标。

$ {\phi }_{1},{\phi }_{2},{\phi }_{3} $$ {\phi }_{4} $正负符号相同,即同为正号或同为负号,则说明接触点在四边形面内,颗粒与四边形单元的接触模式为面接触。

1.3 超二次曲面颗粒的非线性接触力模型

颗粒间及颗粒与结构间的作用力均采用非线性接触力模型[27],其接触力主要由法向力$ {{\boldsymbol{F}}}_{\mathrm{n}} $和切向力$ {{\boldsymbol{F}}}_{\mathrm{t}} $组成。其中,法向力$ {{\boldsymbol{F}}}_{\mathrm{n}} $包括弹性力$ {{\boldsymbol{F}}}_{\mathrm{n}}^{\mathrm{e}} $和黏滞力$ {{\boldsymbol{F}}}_{\mathrm{n}}^{\mathrm{d}} $,可分别表示为

$ {{\boldsymbol{F}}}_{\mathrm{n}}^{\mathrm{e}}=\frac{4}{3}{E}^{*}\sqrt{{r}^{*}}{\left({{\boldsymbol{\delta }}}_{\mathrm{n}}\right)}^{\frac{3}{2}} $ (8)
$ {{\boldsymbol{F}}}_{\mathrm{n}}^{\mathrm{d}}={C}_{\mathrm{n}}{\left({8m}^{*}{E}^{*}\sqrt{{r}^{*}\left|{{\boldsymbol{\delta }}}_{\mathrm{n}}\right|}\right)}^{\frac{1}{2}}\cdot {{{\boldsymbol{v}}}}_{\mathrm{n}} $ (9)

式中:$ E,{C}_{\mathrm{n}},m,r $$ {{{\boldsymbol{v}}}}_{n} $分别为弹性模量、法向阻尼系数、质量、曲率半径和法向相对速度;$ {E}^{*} $$ {r}^{*} $$ {m}^{*} $分别为等效弹性模量、等效曲率半径和等效质量,由$\left(\dfrac{1}{{E}^{\mathrm{*}}}=\dfrac{1-{v}_{i}^{2}}{{E}_{i}}+\dfrac{1-{v}_{j}^{2}}{{E}_{j}},\dfrac{1}{{m}^{\mathrm{*}}}=\dfrac{1}{{m}_{i}}+\dfrac{1}{{m}_{j}}\right.$$\left.\dfrac{1}{{r}^{\mathrm{*}}}=\right. \left. \dfrac{1}{{r}_{i}}+\dfrac{1}{{r}_{j}}\right)$计算得出;$ i\mathrm{和}j $分别表示互相接触的两个单元编号。其中,曲率半径可表示为[27]

$ \begin{split} r=&1\Bigg/\left[\nabla {\boldsymbol{F}}^{\mathrm{T}}{\text{·}} {\nabla }^{2}F{\text{·}} \nabla F-{\left|\nabla F\right|}^{2}\right.\\ &\left.\left(\frac{{\partial }^{2}F}{\partial {\boldsymbol{x}}^{2}}+\frac{{\partial }^{2}F}{\partial {y}^{2}}+\frac{{\partial }^{2}F}{\partial {z}^{2}}\right)\right]\Bigg/2{\left|\nabla F\right|}^{3} \end{split} $ (10)

考虑Mohr-Coulomb摩擦准则,超二次曲面颗粒间的切向力$ {{\boldsymbol{F}}}_{\mathrm{t}} $包括弹性力$ {{\boldsymbol{F}}}_{\mathrm{t}}^{\mathrm{e}} $和黏滞力$ {{\boldsymbol{F}}}_{\mathrm{t}}^{\mathrm{d}} $可分别表示为

$\begin{split} {{\boldsymbol{F}}}_{\mathrm{t}}^{\mathrm{e}}=&{\mu }_{{\rm{s}}}\left|{{\boldsymbol{F}}}_{{\rm{n}}}^{{\rm{e}}}\right|\left\{1-{\left[1-\mathrm{m}\mathrm{i}\mathrm{n}\left({{\boldsymbol{\delta }}}_{\mathrm{t}},{{\boldsymbol{\delta }}}_{\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}\right)\Big/{{\boldsymbol{\delta }}}_{\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}\right]}^{\tfrac{3}{2}}\right\}\cdot\\& {{\boldsymbol{\delta }}}_{\mathrm{t}}\Big/\left|{{\boldsymbol{\delta }}}_{\mathrm{t}}\right| \end{split}$ (11)
$ \begin{split} {{\boldsymbol{F}}}_{\mathrm{t}}^{\mathrm{d}}=&{C}_{\mathrm{t}}\left(6{\mu }_{\mathrm{s}}{{\rm{m}}}^{*}\left|{{\boldsymbol{F}}}_{{\rm{n}}}^{{\rm{e}}}\right|\cdot \right. \\& \left. \sqrt{1-\mathrm{m}\mathrm{i}\mathrm{n}\left({{\boldsymbol{\delta }}}_{\mathrm{t}},{{\boldsymbol{\delta }}}_{\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}\right)/{{\boldsymbol{\delta }}}_{{\rm{t}},\mathrm{m}\mathrm{a}\mathrm{x}}}\right)^{\tfrac{1}{2}}\cdot {{{\boldsymbol{v}}}}_{\mathrm{t}} \end{split} $ (12)

式中:${\mu }_{\mathrm{s}},{C}_{\mathrm{t}}$${{\boldsymbol{v}}}_{\mathrm{t}}$分别为滑动摩擦系数、切向阻尼系数和切向相对速度;${{\boldsymbol{\delta }}}_{\mathrm{t}}={{\boldsymbol{\delta }}}_{\mathrm{t}}+{\boldsymbol{v}}_{\mathrm{t}}{\varDelta }_{{\rm{t}}}$为切向相对位移;$ {\varDelta }_{{\rm{t}}} $为DEM时间步长;$ {{\boldsymbol{\delta }}}_{\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}} $为最大切向重叠量,可表示为$ {{\boldsymbol{\delta }}}_{\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}={\mu }_{\mathrm{s}}\left(2- \nu\right)/2\left(1- \nu\right){{\boldsymbol{\delta }}}_{\mathrm{n}} $$ \nu $为泊松比。

1.4 超二次曲面颗粒的信息更新

超二次曲面颗粒的运动方程是基于牛顿第二定律,每个颗粒在当前时间步的平动速度和旋转速度可表示为

$ m\frac{\mathrm{d}{{\boldsymbol{v}}}_{t}}{\mathrm{d}t}=\sum _{i=1}^{N}\left({{\boldsymbol{F}}}_{\mathrm{n},i}+{{\boldsymbol{F}}}_{\mathrm{t},i}\right)+m{{\boldsymbol{g}}} $ (13)
$ {{\boldsymbol{J}}}\frac{{\mathrm{d}{{\boldsymbol{\omega}} }}_{t}}{\mathrm{d}t}=\sum _{i=1}^{N}\left({{{\boldsymbol{M}}}}_{\mathrm{n},i}+{{{\boldsymbol{M}}}}_{\mathrm{t},i}+{{{\boldsymbol{M}}}}_{\mathrm{r},i}\right) $ (14)

式中:$ {{\boldsymbol{v}}} $$ {{\boldsymbol{\omega}} } $${\boldsymbol{m}} $$ {{\boldsymbol{J}}} $分别为颗粒的平动速度、转动速度、质量和转动惯量;${{\boldsymbol{F}}}_{{\mathrm{n}},i}$${{\boldsymbol{F}}}_{{\mathrm{t}},i}$分别为第$ i $个接触对的法向力和切向力;$ {{{\boldsymbol{M}}}}_{\mathrm{n},i} $$ {{{\boldsymbol{M}}}}_{\mathrm{t},i} $$ {{{\boldsymbol{M}}}}_{\mathrm{r},i} $分别为第$ i $个接触对的法向力和切向力产生的力矩以及滚动摩擦力矩;$ N $为总接触对数;$ {{\boldsymbol{g}}} $为重力加速度。

超二次曲面颗粒在旋转时通常会引起惯性张量的改变。因此,引入每个颗粒的局部坐标系来反映颗粒的转动,同时采用四元数方法计算全局坐标和局部坐标间的转化关系,可表示为

$ {{{\boldsymbol{e}}}}_{\mathrm{l}}^{\mathrm{T}}={{\boldsymbol{R}}}{{{\boldsymbol{e}}}}_{\mathrm{g}}^{\mathrm{T}} $ (15)

式中,${\boldsymbol{e}}_{\mathrm{l}}^{\mathrm{T}}$${\boldsymbol{e}}_{\mathrm{g}}^{\mathrm{T}}$分别为局部坐标系下的方向矢量和全局坐标系下的方向矢量。

1.5 有限元方法中六面体实体单元的信息更新

颗粒与结构间的接触力模型与颗粒间的接触力模型相同。它根据颗粒与结构的重叠量计算相应的接触力,并将接触力转化为有限元分析中的等效节点力。根据虚功定理,原始接触力的虚功之和等于所有等效节点力的虚功之和。对于一个包含8个节点的实体单元,其等效节点力可以表示为

$ {\left\{{{\boldsymbol{F}'}}\right\}}={\left({N}\right)}_{j}\left\{{{\boldsymbol{f}}}\right\} $ (16)

式中:${\left\{{{\boldsymbol{F}'}}\right\}}$是等效节点力;${\left({N}\right)}_{j}$为在单元$j$的形函数;$\left\{{{\boldsymbol{f}}}\right\}$为接触力。

有限元方法是一种有效的数值计算方法,通过在每个单元内对位移进行逼近,将连续体转化为有限自由度的系统,用以描述结构体的运动行为。在本文中,通过采用8节点六面体实体单元的线弹性有限元方法建立工程结构的动力模型。其中,8节点六面体实体单元的单元刚度矩阵可表示为

$ {{\boldsymbol{K}}}={\int }_{-1}^{1}{\int }_{-1}^{1}{\int }_{-1}^{1}{{{\boldsymbol{B}}}}^{\mathrm{T}}{{\boldsymbol{D}}}{{\boldsymbol{B}}}\left|{{\boldsymbol{J}}}\right|\mathrm{d}\xi \mathrm{d}\zeta \mathrm{d}\eta $ (17)

式中:$ {{\boldsymbol{K}}} $为六面体单元的单元刚度矩阵;$ {{\boldsymbol{B}}} $为单元的几何矩阵;$ {{\boldsymbol{D}}} $为单元的弹性矩阵;$ {{\boldsymbol{J}}} $为雅可比矩阵;$\xi ,\zeta$$ \eta $分别为插值点的自然坐标。

为了提高有限元分析的计算效率,本研究中采用了一种基于EBE (element-by-element)的隐式积分方法,用于求解有限元动力学方程。有限元动力学方程可以表示为

$ {{\boldsymbol{M}}}{\ddot{{{\boldsymbol{u}}}}}_{t}+{{\boldsymbol{C}}}{\dot{{{\boldsymbol{u}}}}}_{t}+{{\boldsymbol{K}}}{{{\boldsymbol{u}}}}_{t}={{\boldsymbol{F}}}_{t} $ (18)

式中:$ {\ddot{{{\boldsymbol{u}}}}}_{t} $$ {\dot{{{\boldsymbol{u}}}}}_{t} $$ {{{\boldsymbol{u}}}}_{t} $分别表示$ t $时刻结构的加速度、速度和位移;$ {{\boldsymbol{M}}} $$ {{\boldsymbol{K}}} $分别为结构的集中质量矩阵和刚度矩阵;$ {{\boldsymbol{F}}}_{t} $$ t $时刻受到的外部载荷;$ {{\boldsymbol{C}}} $为阻尼矩阵。

2 DEM-FEM耦合算法的理论验证

为验证超二次曲面DEM-FEM耦合算法的可靠性,这里模拟了薄板受颗粒重力作用时的位移随长度的变化,并将数值结果与基于梁挠度方程的解析结果以及有限元ANSYS软件的数值结果进行对比验证。采用超二次曲面方程构造30000个不同形态的颗粒,包括10000个立方体$\left(a=b=c=\right. \left.{3\;\mathrm{m}\mathrm{m}\mathrm{和}n}_{1}={n}_{2}=6\right)$、10000个圆柱体($a=b=3\;\mathrm{m}\mathrm{m}\mathrm{,} {c=4\;\mathrm{m}\mathrm{m}\mathrm{,}n}_{1}=8\mathrm{和}{n}_{2}=2$)和10000个球体($a=b= c= 3\;\mathrm{m}\mathrm{m}$$ {n}_{1}={n}_{2}=2 $)。薄板上方所有颗粒以任意位置生成并在平板上实现稳定堆积,如图3所示。


图 3 球体和非球体颗粒在薄板上的稳定堆积 Fig. 3 Stable packing of the spherical and non-spherical particles on the thin plane

以左侧薄板中心为坐标原点,对薄板在$ x=0 $处进行位移约束,即${u}_{x}=0\mathrm{,}{u}_{y}=0\mathrm{和}{u}_{z}=0$。DEM-FEM耦合计算参数列于表1中。在有限元ANSYS软件中,将所有颗粒的重力等效为薄板上的均布面力,并且对薄板在$ x=0 $处所有自由度进行约束。考虑一端受全自由度约束且受均布载荷作用下梁的挠度方程可表示为


表 1 DEM-FEM耦合模拟的主要参数 Table 1 Main parameters of DEM-FEM coupled simulation
$ y=-\frac{q{\boldsymbol{x}}^{2}}{24EI}\left({\boldsymbol{x}}^{2}-4lx+6{l}^{2}\right) $ (19)

式中:$ q $为作用在梁上的均布载荷,$ q=\dfrac{{N}_{\mathrm{l}}{m}_{\mathrm{l}}g}{S} $$ {N}_{\mathrm{l}} $$ {m}_{\mathrm{l}} $分别为颗粒数目和颗粒质量,$ S $为薄板的横截面积;$ I $为薄板的转动惯量;$ l $为薄板的长度。

当颗粒材料无相对运动时,对薄板的垂直位移${d}_{z}$进行统计,并将DEM-FEM耦合算法得出的数值结果与ANSYS软件得到的结果进行对比,如图4所示。结果显示,DEM-FEM耦合算法得出的位移大小和分布规律与ANSYS软件得到的数值结果基本一致。同时,将计算得到的薄板垂直位移随长度变化的数值结果与式(19)的理论结果及ANSYS软件得到的数值结果进行对比,如图5所示。薄板的垂直位移随着长度的增加而减小,且下降速率逐渐增加。ANSYS软件的计算结果与梁的解析结果比较吻合,而DEM-FEM耦合算法得到的数值结果与解析结果存在微小差异。这些微小差异主要是由于薄板上的颗粒位置具有随机性,这与解析结果中的均布荷载存在一定区别。


图 4 由DEM-FEM耦合算法和ANSYS软件得到的薄板位移云图对比 Fig. 4 Comparison of the displacement patterns of a thin plate obtained by the DEM-FEM coupled algorithm and the ANSYS software

图 5 薄板竖直方向位移随薄板长度的变化关系 Fig. 5 Relationship between the vertical displacement of the plate and the length
3 基于DEM-FEM耦合算法的滚筒球磨机模拟

采用基于超二次曲面方程的非规则颗粒DEM-FEM耦合算法模拟滚筒球磨机中不同形态颗粒的运动行为,颗粒形状包括球体、圆柱体和立方体颗粒。滚筒球磨机的内壁半径${r}_{\mathrm{m}}=0.25\;\mathrm{m}$,筒壁厚度${w}_{\mathrm{m}}=0.01\;\mathrm{m}$,提升条长度${s}_{\mathrm{l}}=0.26\;\mathrm{m}$,高度${h}_{\mathrm{l}}=0.26\;\mathrm{m}$,整体宽度${d}_{\mathrm{m}}=0.25\;\mathrm{m}$,如图6所示。滚筒球磨机模型共有2880个实体单元和4752个节点,同时考虑对提升条进行局部网格加密[28]。颗粒密度为$ 7\;900\;\mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $,时间步长为2.0$ {\times 10}^{-6}\;\mathrm{s} $,其他DEM-FEM模拟参数列于表1中。在有限元计算球磨机位移及应力分布时,球磨机外侧筒壁的自由度将被全部约束。


图 6 滚筒球磨机的有限元模型 Fig. 6 Finite element modeling of the tumbling mill

采用超二次曲面方程分别构造2000个球体($ {n}_{1}={n}_{2}=2 $)、圆柱体($ {n}_{1}=8,{n}_{2}=2 $)和立方体($ {n}_{1}= {n}_{2}=6 $)。所有颗粒具有相等质量,半轴长比例均满足$ a:b:c=1:1:1 $,等效半径为5.8 mm。在初始时刻,所有颗粒具有随机位置和空间角度。在重力作用下,颗粒自由下落至滚筒内并形成稳定堆积。当所有颗粒无相对运动时,滚筒以转速为30 r/min开始旋转。图7为滚筒球磨机中球体、圆柱体和立方体颗粒在不同时刻的运动状态。4幅插图分别对应颗粒材料在滚筒球磨机中生成、下落稳定、球磨机旋转$ {90}\text{°} $和球磨机旋转$ {180}\text{°} $时的运动状态。


图 7 滚筒球磨机中不同形态颗粒材料的运动状态 Fig. 7 Motion of the granular materials with different morphologies in the tumbling mill

图8为滚筒球磨机中球体、圆柱体和立方体颗粒材料的速度分布。其中,V为颗粒在3个主轴方向的合速度。相较于球体,圆柱体和立方体在颗粒提升过程中形成了更为稳定的颗粒运动。其中,立方体颗粒被提升的位置最高。这是因为这些颗粒具有更大的相对速度,这导致颗粒间碰撞更加剧烈。因此,立方体颗粒材料具有最好的研磨效果,而圆柱体颗粒材料比球体颗粒材料具有更好的研磨效果。


图 8 不同形态颗粒在滚筒球磨机中的速度分布图 Fig. 8 Velocity distribution of particles with different morphology in the tumbling mill

随着滚筒球磨机的旋转,颗粒在重力作用下不断与球磨机壁面发生接触和碰撞,这种荷载作用使得结构在接触点处产生局部变形,并进一步产生应力–应变响应。图9为滚筒球磨机在颗粒碰撞下的局部von-mises应力云图。结果表明:滚筒球磨机中最大应力出现在提升条上。提升条的作用是提供颗粒克服自身重力的支持力,并不断提升颗粒的高度。因此,颗粒与球磨机间的剧烈碰撞主要发生在提升条上,这与之前文献[28]的研究结果是一致的。


图 9 滚筒球磨机受载荷部分的von-mises应力云图 Fig. 9 Von-mises stress pattern of the loaded part in the tumbling mill

当球体颗粒材料填充进滚筒球磨机时,提升条位移${d}_{\mathrm{l}}$随球磨机旋转角度φ的变化如图10所示。在球磨机的运动过程中,提升条位移峰值的波动较小,约为$ 8.2\times {10}^{-7}\;\mathrm{m} $。这些峰值位移分别出现在$ 4{0}\text{°}~11{0}\text{°} $$ 32{0}\text{°}~47{0}\text{°} $以及$ {670}\text{°}~{700}\text{°} $。这些角度范围分别对应了球磨机旋转的第一周中颗粒最初抬升、第一周末尾的颗粒自由下落和第二周中颗粒再次抬升,以及第二周末尾颗粒再次自由下落的阶段。从位移变化趋势可以看出:起初位移迅速上升,达到峰值后进入波动阶段;随后,位移峰值不断减小,直到约$ 11{0}\text{°} $时进入无负载状态,位移迅速减小并稳定下来;球磨机旋转至约$ 30{0}\text{°} $时再次进入负载状态,位移开始以较小周期和较大振幅波动,直到约$ 47{0}\text{°} $时再次进入第二周的无负载状态;此外,球磨机在约$ {670}\text{°} $时再次进入负载状态。值得注意的是,这一阶段的变化特性与在$ 30{0}\text{°}~36{0}\text{°} $区间的特性基本相似,这也说明当球磨机进入稳定状态后,其有限元位移特征基本相同。


图 10 球体颗粒填充下提升条位移随旋转角度的变化关系 Fig. 10 Relationship between the displacement of lifting bar and rotation angle under spherical granular filling

当圆柱体颗粒材料填充进滚筒球磨机时,提升条位移${d}_{\mathrm{l}}$随球磨机旋转角度φ的变化如图11所示。第一周颗粒抬升阶段的位移峰值约为$ 6.3\times {10}^{-7}\;\mathrm{m} $,而第二周位移峰值增加到$ 1.0\times {10}^{-6}\;\mathrm{m} $。与球体颗粒材料填充时相比,圆柱体颗粒材料填充时的位移大于初始提升阶段的位移。从图8(b)中可以看出,圆柱体颗粒材料随着球磨机的旋转而逐渐被提升且自由下落至颗粒床表面。这加剧了颗粒与颗粒、颗粒与球磨机间的碰撞,使下落后的颗粒获得更大的动能,从而使稳定阶段的位移更大。与球体颗粒材料相似,圆柱体颗粒材料的位移也相应出现了3个阶段,对应的角度范围分别为$ {0}\text{°} $${120}\text{°} $$ {305}\text{°} $${490}\text{°} $以及$ {670}\text{°} $$72{0}\text{°} $。角度范围不同的原因主要在于圆柱体颗粒比球体颗粒具有更加不规则的表面,进而导致颗粒间产生更强的互锁,并使得颗粒在提升条作用下被提升至更高位置。从位移的变化趋势上可以看出,两次颗粒提升阶段的位移峰值都出现在最后一个阶段,分别约为$ {100}\text{°} $$ {450}\text{°} $,并且波动的振幅小于球体颗粒填充时得到的位移幅值。这是由于圆柱体颗粒的运动更加稳定,而球体颗粒更容易发生滑动和滚动,导致提升条受到的载荷变化更加剧烈。


图 11 圆柱体颗粒填充下提升条位移随旋转角度的变化关系 Fig. 11 Relationship between the displacement of lifting bar and rotation angle under cylindrical granular filling

当立方体颗粒材料填充进滚筒球磨机时,提升条位移${d}_{\mathrm{l}}$随球磨机旋转角度φ的变化如图12所示。提升条产生位移的角度范围为$ {0}\text{°} $$ {130}\text{°} $$ {330}\text{°} $${500}\text{°} $$ {675}\text{°} $${720}\text{°} $区间。与圆柱体颗粒材料填充时提升条的位移变化类似,提升条在稳定阶段的位移峰值大于初始抬升阶段的位移峰值,且峰值约为$ 7.1\times {10}^{-7}\;\mathrm{m} $$ 1.2\times {10}^{-6}\;\mathrm{m} $,两处峰值均大于圆柱体颗粒材料填充时提升条的位移峰值。在旋转的第一周末尾,即角度范围在$ {300}\text{°}~{360}\text{°} $,立方体颗粒材料对提升条产生的位移明显小于圆柱体和球体颗粒对提升条产生的位移。这是因为立方体颗粒比圆柱体和球体颗粒具有更尖锐的表面和更强的互锁效应,这使得立方体颗粒可以更好地借助提升条向上运动,同时,在旋转末尾阶段产生颗粒的滞后下落现象。


图 12 立方体颗粒填充下提升条位移随旋转角度的变化关系 Fig. 12 Relationship between the displacement of lifting bar and rotation angle under cubic granular filling
4 结束语

采用基于连续函数包络的超二次曲面方程构造非球形颗粒材料,在传统球体颗粒的离散元–有限元耦合算法基础上,发展了非球体颗粒材料的离散元–有限元耦合算法。该算法采用超二次曲面颗粒与8节点实体单元间的界面搜索方法计算接触力,同时通过形函数将接触力插值至节点上。随后,基于等效节点力求解有限元方法中的隐式积分动力方程,并更新超二次曲面颗粒及结构的位置信息。为验证基于超二次曲面的离散元–有限元耦合算法的有效性,对由颗粒重力引起的薄板变形进行统计,其薄板位移的数值结果与解析结果和有限元ANSYS软件的数值结果均有较好的吻合。由此表明:基于超二次曲面的离散元–有限元耦合算法能够合理描述非球体颗粒与结构之间的相互作用,同时也能准确反映结构的位移变化。在此基础上,采用超二次曲面方程构造了球体、圆柱体和立方体颗粒,并对滚筒球磨机中不同形态颗粒材料的流动行为、速度分布、结构位移变化以及局部应力分布进行数值分析,该研究结果可为滚筒球磨机的内部结构设计提供重要参考。

参考文献
[1]
SHAHBAZI B, JAFARI M, PARIAN M, et al. Study on the impacts of media shapes on the performance of tumbling mill—A review[J]. Minerals Engineering, 2020, 157: 106490. DOI:10.1016/j.mineng.2020.106490
[2]
NKWANYANA S, LOVEDAY B. Addition of pebbles to a ball-mill to improve grinding efficiency–Part 2[J]. Minerals Engineering, 2018, 128: 115-122. DOI:10.1016/j.mineng.2018.08.024
[3]
SHI F. Comparison of grinding media—Cylpebs versus balls[J]. Minerals Engineering, 2004, 17(11/12): 1259-1268.
[4]
QIAN H Y, KONG Q G, ZHANG B L. The effects of grinding media shapes on the grinding kinetics of cement clinker in ball mill[J]. Powder Technology, 2013, 235: 422-425. DOI:10.1016/j.powtec.2012.10.057
[5]
SIMBA K P, MOYS M H. Effects of mixtures of grinding media of different shapes on milling kinetics[J]. Minerals Engineering, 2014, 61: 40-46. DOI:10.1016/j.mineng.2014.03.006
[6]
CLEARY P W. Effect of rock shape representation in DEM on flow and energy utilisation in a pilot SAG mill[J]. Computational Particle Mechanics, 2019, 6(3): 461-477. DOI:10.1007/s40571-019-00226-3
[7]
VERMEULEN L A, HOWAT D D. A sampling procedure validated[J]. Journal of the South African Institute of Mining and Metallurgy, 1989, 89(12): 365-370.
[8]
LU G, THIRD J R, MÜLLER C R. Discrete element models for non-spherical particle systems: from theoretical developments to applications[J]. Chemical Engineering Science, 2015, 127: 425-465. DOI:10.1016/j.ces.2014.11.050
[9]
陈泉, 杨晖, 李然, 等. 转筒颗粒流崩塌前后的颗粒温度和体积分数测量[J]. 中国科学:物理学 力学 天文学, 2019, 49(6): 067001.
[10]
童红政, 袁岑颉, 薛晓垒, 等. 中速磨煤机构件对气相场影响的数值模拟[J]. 能源研究与信息, 2023, 39(3): 166-4759. DOI:10.13259/j.cnki.eri.2023.03.005
[11]
常晓林, 马刚, 周伟, 等. 颗粒形状及粒间摩擦角对堆石体宏观力学行为的影响[J]. 岩土工程学报, 2012, 34(4): 646-653.
[12]
CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47-65.
[13]
郭宇, 凡凤仙, 白鹏博, 等. 颗粒物质在竖直振动U形管中迁移的离散元方法模拟[J]. 上海理工大学学报, 2019, 41(5): 409-416. DOI:10.13255/j.cnki.jusst.2019.05.001
[14]
许文祥, 孙洪广, 陈文, 等. 软物质系颗粒材料组成、微结构与传输性能之间关联建模综述[J]. 物理学报, 2016, 65(17): 178101.
[15]
孙其诚, 王光谦. 颗粒流动力学及其离散模型评述[J]. 力学进展, 2008, 38(1): 87-100.
[16]
GAN J Q, ZHOU Z Y, YU A B. Interparticle force analysis on the packing of fine ellipsoids[J]. Powder Technology, 2017, 320: 610-624. DOI:10.1016/j.powtec.2017.07.064
[17]
WANG S Q, FAN Y N, JI S Y. Interaction between super-quadric particles and triangular elements and its application to hopper discharge[J]. Powder Technology, 2018, 339: 534-549. DOI:10.1016/j.powtec.2018.08.026
[18]
刘璐, 季顺迎. 基于扩展多面体包络函数的快速接触搜索算法[J]. 中国科学:物理学 力学 天文学, 2019, 49(6): 064601.
[19]
ZHONG W Q, YU A B, LIU X J, et al. DEM/CFD-DEM modelling of non-spherical particulate systems: theoretical developments and applications[J]. Powder Technology. 2016, 302: 108–152.
[20]
崔泽群, 陈友川, 赵永志, 等. 基于超二次曲面的非球形离散单元模型研究[J]. 计算力学学报, 2013, 30(6): 854-859. DOI:10.7511/jslx201306017
[21]
CLEARY P W, SINNOTT M D, MORRISON R D, et al. Analysis of cone crusher performance with changes in material properties and operating conditions using DEM[J]. Minerals Engineering, 2017, 100: 49-70. DOI:10.1016/j.mineng.2016.10.005
[22]
WANG S Q, ZHOU Z Y, JI S Y. Radial segregation of a gaussian-dispersed mixture of superquadric particles in a horizontal rotating drum[J]. Powder Technology, 2021, 394: 813-824. DOI:10.1016/j.powtec.2021.09.012
[23]
ZHENG Z M, ZANG M Y, CHEN S H, et al. A GPU-based DEM-FEM computational framework for tire-sand interaction simulations[J]. Computers & Structures, 2018, 209: 74-92.
[24]
BARR A H. Superquadrics and angle-preserving transformations[J]. IEEE Computer Graphics and Applications, 1981, 1(1): 11-23. DOI:10.1109/MCG.1981.1673799
[25]
ZHAO Y Z, XU L, UMBANHOWAR P B, et al. Discrete element simulation of cylindrical particles using super-ellipsoids[J]. Particuology, 2019, 46: 55-66. DOI:10.1016/j.partic.2018.04.007
[26]
PODLOZHNYUK A, PIRKER S, KLOSS C. Efficient implementation of superquadric particles in Discrete Element Method within an open-source framework[J]. Computational Particle Mechanics, 2016, 4(1): 101-118.
[27]
王嗣强, 季顺迎. 考虑等效曲率的超二次曲面单元非线性接触模型[J]. 力学学报, 2018, 50(5): 1081-1092. DOI:10.6052/0459-1879-18-103
[28]
JONSÉN P, PÅLSSON B I, TANO K, et al. Prediction of mill structure behaviour in a tumbling mill[J]. Minerals Engineering, 2011, 24(3/4): 236-244.