2. 精密超精密制造装备及控制北京市重点实验室, 北京 100084
2. Beijing Key Lab of Precision/Ultra-precision Manufacturing Equipments and Control, Beijing 100084, China
蛋白质工程是当今科学研究的重点。目前，获取蛋白质的方法主要有两大类： 一类是通过化学方法人工合成，另一类则是直接从生物组织中提取。其中从生物组织中提取出的蛋白质是人工合成的范例与基础。从生物组织提取蛋白质的第一个工序是将生物组织粉碎成微米级的粉末。生物样品粉末的粒度大小直接决定蛋白质的提取率。在生物组织粉碎中，振动球磨机构应用最为广泛。它的工作原理是利用研磨罐的振动激励罐内研磨球(多数为钢球或石球)高速、高频碰撞产生的挤压力、剪切力、摩擦力等将磨料粉碎。
振动球磨机构结构简单，但研磨罐内的粉碎过程却十分复杂。为了解球磨机构粉碎工作机理和研磨行为，很多学者对研磨罐内的研磨球进行了运动学和动力学研究。离散元法(discrete element method,DEM)常被用于模拟球磨机构中的研磨过程，通过该模型可以直观地了解磨球运动，分析研磨条件对研磨过程的影响[1, 2, 3, 4, 5, 6, 7]。然而，因为离散元法的每一步迭代都需要计算作用在磨球上的力并更新磨球的速度和位置，所以计算量大、耗时长。因而，有学者提出运用事件驱动法研究研磨过程[8, 9, 10, 11, 12]。Gavrilov等[8]采用事件驱动法建立了研磨介质在研磨罐内的动力学模型。经过对立式振动球磨仪的动力学仿真，分析了振动频率及振动幅度与耗散能量的关系，计算了磨球碰撞次数以及耗散能量，结合粉体的随机破碎模型模拟了研磨过程中粉体粒径的变化。为了简化计算，该模型假定： 碰撞均为正面碰撞，且忽略磨球自转。但是，实际运动中磨球带有角速度，且其角速度对磨球碰撞具有较大的影响，另外，已有研究表明，振动球磨机构中50%以上的碰撞角度为80°~90°，故大多数属于斜碰[1, 2]。Reichardt等[11]用事件驱动法对滚筒式球磨仪进行动力学仿真，并分析了法向和切向恢复系数对磨球速度的影响。为了验证动力学模型的正确性，文[8, 11]都仿真了理想气体模型中颗粒的速度分布，得到的结果较接近于Maxwell速度分布。因为事件驱动法只考虑颗粒系统中的颗粒与颗粒碰撞发生的时刻，不考虑在2个碰撞事件之间所发生的力学行为，所以仿真速度快、效率高。该方法简单计算量小，故更适合于稀疏颗粒系统的仿真与建模。
本文针对清华大学机器人及其自动化实验室研制的蛋白质自动粗分离仪器(Tsinghua protein extraction device,THPED)中的立式振动球磨机构进行研究，分析研磨参数对研磨过程特征参量的影响。首先，在考虑磨球自转的情况下，采用事件驱动法仿真分析单个磨球在研磨罐内的动力学特性，得出磨球与研磨罐壁之间的碰撞频率(即单位时间内磨球与研磨罐壁之间的碰撞次数)和碰撞能量两种特征参量。然后分析总结这两种特征参量随振动频率、振动幅度以及研磨罐直径变化的规律，为研磨机构设计和参数选择提供理论依据。
1 研磨机构及仿真方法 1.1 研磨机构简介清华大学机器人及其自动化实验室研制的蛋白质自动粗分离仪器将研磨与离心机构集成在一起，能从生物组织中自动提取出蛋白质。其研磨机构见图1，主要由研磨罐、从动杆、基座、直线轴承及驱动凸轮5部分组成。其中，研磨罐的磨腔由两部分组成： 中间是一个直径为Dv、 高度为H的圆柱形，两端为半球形。
|
| 图1 THPED的研磨机构简图 |
驱动凸轮是一个由电机驱动的圆柱凸轮，凸轮上的曲槽是一个双周期的正弦曲线。当凸轮旋转时驱动从动杆连同研磨罐一起做直线往复运动，研磨罐内的钢球与研磨罐发生碰撞，从而使样品得到充分地粉碎。对于不同的应用场合要求，THPED只要改变电机的转速就能方便地调整振动频率。
1.2 仿真方法事件驱动法常用于颗粒系统的建模与仿真[13, 14, 15]。该方法采用非弹性硬球模型(inelastic hard sphere model)，即所有碰撞都被认为是瞬间发生的，颗粒可被看作是一个硬球，碰撞时不引起颗粒本身的弹性变形，碰撞过程中的能量损失可通过恢复系数表示。事件驱动法实际上是一个寻找颗粒系统中发生碰撞的时刻以及处理碰撞响应的过程，由于粉体的质量和尺寸远小于磨球的质量和尺寸，故可忽略粉体对磨球运动和碰撞的影响，因此本文的仿真模型只研究磨球和研磨罐2个对象。采用事件驱动法对单个磨球在研磨罐内的动力学仿真基本流程见图2，包括3个步骤：
|
| 图2 事件驱动法仿真流程 |
步骤1 预计磨球与研磨罐的碰撞时刻。假设磨球在研磨罐内的运动不受自身角速度以及空气阻力的影响，即磨球只受到重力加速度的作用，那么磨球在研磨罐内做斜抛运动。磨球与研磨罐内壁的碰撞分3种可能，分别是磨球与中间圆柱面碰撞、磨球与上半球面碰撞以及磨球与下半球面碰撞。根据研磨罐和磨球的运动学方程以及研磨罐内壁的几何形状分别建立出3个求解碰撞时刻的方程，3个方程最小正解作为磨球与研磨罐内壁预计碰撞时刻。
步骤2 计算碰撞后磨球的速度。由于非弹性硬球模型认为碰撞是瞬间发生的，所以碰撞发生后，磨球碰撞前速度立即被碰撞后速度代替。磨球碰撞前后的速度关系为
u′=Ψ(u,r).
步骤3 更新磨球位置和速度。根据前2步计算结果更新磨球位置和速度，准备进入下一个循环。
2 动力学模型 2.1 碰撞检测碰撞检测目的是确定磨球与研磨罐发生碰撞的时刻，首先建立坐标系如图3所示。其中，坐标系原点与初始化位置的研磨罐中心重合，Z轴方向与研磨罐上下运动方向重合，坐标系遵守右手法则。
|
| 图3 磨球与研磨罐碰撞 |
研磨罐中心和磨球的运动方程分别为
\[\left\{ \begin{matrix}
X=Y=0, \\
Z=A\sin \left( 2\pi ft \right) \\
\end{matrix} \right.,\]
(1)
\[{{r}_{b}}={{r}_{0}}+{{v}_{b}}\Delta t+\frac{1}{2}g\Delta {{t}^{2}}.\]
(2)
如节1.2中所述，磨球与研磨罐的碰撞分3种可能。在磨球与中间圆柱面碰撞的情况下，只需要考虑磨球在水平面OXY中的位置。若球与研磨罐的圆柱面发生碰撞，则球心到圆柱面轴线距离应为： R=Rv－Rb，Rv和Rb分别为研磨腔和磨球的半径。再结合式、得到碰撞条件为
\[\left\{ \begin{matrix}
r_{Xb}^{2}+r_{Yb}^{2}={{R}^{2}}, \\
|{{r}_{Zb}}-Asin\left( 2\pi ft \right)| \\
\end{matrix} \right.\le \frac{H}{2}.\]
(3)
rXb=rXb0+vXbΔt，rYb=rYb0+vYbΔt，rZb=rZb0+vZbΔt－$\frac{1}{2}$gΔt2.
设式(3)的解为tcyl。若tcyl不大于t0，则令tcyl=+∞。
磨球与研磨罐下半球面碰撞时，磨球心到曲面球心的距离应为R=Rv－Rb，并且磨球心位置应低于下半球面的球心位置，从而得到碰撞条件为
\[\left\{ \begin{matrix}
r_{Xb}^{2}+r_{Yb}^{2}+{{\left( {{r}_{Zb}}-Asin\left( 2\pi ft \right)+\frac{H}{2} \right)}^{2}}={{R}^{2}}, \\
{{r}_{Zb}}<-\frac{H}{2}+Asin\left( 2\pi ft \right). \\
\end{matrix} \right.\]
(4)
类似地，磨球与研磨罐的上半球面碰撞的碰撞时刻是式(5)的解。
\[\left\{ \begin{matrix}
r_{Xb}^{2}+r_{Yb}^{2}+{{\left( {{r}_{Zb}}-Asin\left( 2\pi ft \right)-\frac{H}{2} \right)}^{2}}={{R}^{2}}, \\
{{r}_{Zb}}>-\frac{H}{2}+Asin\left( 2\pi ft \right). \\
\end{matrix} \right.\]
(5)
由于式(4)、 (5)的结构复杂，本文拟用二分法求解。
Z轴方向上，磨球的运动范围总在磨腔的下端点和上端点之间。
磨球达到下端点的下限位置所需要的时间满足：
\[{{r}_{Zb0}}+{{v}_{Zb}}\Delta t-\frac{1}{2}g\Delta {{t}^{2}}=-\left( A+\frac{H}{2}+R \right).\]
\[\frac{{{v}_{Zb}}+\sqrt{{{v}_{Zb}}^{2}+2g\left( {{r}_{Zb0}}+A+\frac{H}{2}+R \right)}}{g}+{{t}_{0}}.\]
(6)
类似地，磨球达到上端点的上限位置所需要的时间为
\[\frac{{{v}_{Zb}}-\sqrt{{{v}_{Zb}}^{2}+2g\left( {{r}_{Zb0}}-A-\frac{H}{2}-R \right)}}{g}+{{t}_{0}}.\]
(7)
碰撞检测时，首先解式(3)得到磨球与研磨罐的圆柱面的碰撞时刻tcyl，然后在(t0,min(tcyl,t1)］区域上用二分法求解式(4)，得到磨球与研磨罐的下半球面的碰撞时刻tlow，最后在(t0,min(tcyl,t2,tlow)］上用二分法求解式(5)得到磨球与研磨罐上半球面的碰撞时刻tup，取3个碰撞时刻中的最小值作为磨球与研磨罐的预计碰撞时刻。
2.2 碰撞响应本文采用非弹性硬球模型作为碰撞模型，磨球与研磨罐接触认为是点接触。因而忽略外力对磨球产生的扭矩。
磨球碰撞后的线速度和角速度分别为[9]
\[v_{b}^{'}={{v}_{b}}-\left( +{{\varepsilon }^{n}} \right)\delta _{bw}^{n}+\frac{2}{7}\left( {{\varepsilon }^{t}}-1 \right)\delta _{bw}^{t}\]
(8)
\[\omega _{b}^{'}={{\omega }_{b}}-\frac{5\left( {{\varepsilon }^{t}}-1 \right)}{7{{R}_{b}}}{{n}_{bw}}\delta _{bw}^{t}.\]
(9)
对于一个振动球磨机构来说，表征其研磨能力的主要特征参量有碰撞频率、碰撞能量、碰撞速度以及碰撞角度等。文[2, 6, 16]表明碰撞频率和碰撞能量是对研磨过程影响最大的参数，碰撞能量越大，产品粒度越小，破碎率越高，因而研磨时间越短。磨球与研磨罐壁碰撞的碰撞能量可如下计算：
\[{{E}_{col}}=\frac{1}{2}{{m}_{b}}v_{bw}^{2}.\]
(10)
为了说明研磨参数如振动频率、振动幅度以及研磨罐直径对研磨特征参量如碰撞频率、碰撞能量的影响，本文进行了一系列仿真实验。对于每一个不同的研磨参数分别做了5组仿真实验，取平均值作为最后结果。仿真时，磨球的初始化位置和速度由程序随机给出，线速度在0~0.1 m/s范围内，角速度在0~0.1 rad/s范围内。仿真过程中用到的参数如下： φ8 mm钢球一个，其材料的密度ρ=7 850 kg/m3，仿真时间为60 s，f=10~25 Hz，A=15~30 mm，法向和切向恢复系数为0.9。对于不同的研磨罐，仿真过程中磨腔体积保持不变(均为5 ml)，只改变其直径大小。研磨罐尺寸如表1所示。
当f=10 Hz、 A=25 mm时，研磨罐直径变化对研磨过程的碰撞频率和碰撞能量的影响如图4所示。随着研磨罐的直径增大，碰撞频率和单次碰撞能量(总碰撞能量与总碰撞次数之比)均减小。因此，设计研磨罐尺寸时，在保证研磨罐体积的条件下，应尽量减小研磨罐的直径。
|
| 图4 不同研磨罐直径下碰撞频率和单次碰撞能量 |
在f=10 Hz、 DV=12 mm时，碰撞频率和单次碰撞能量均与振动幅度成线性关系，如图5所示。
|
| 图5 不同振动幅度下碰撞频率和单次碰撞能量 |
在A=25 mm、 DV=12 mm的条件下，得到碰撞频率随振动频率的变化见图6a，可见二者成线性关系。用二次多项式对实验结果进行拟合得到单次碰撞能量与振动频率之间呈平方关系(见图6b)，这表明单次碰撞能量对振动频率更为敏感。
|
| 图6 不同振动频率下碰撞频率和单次碰撞能量 |
为了验证上述仿真分析结果，本文利用图1所示的THPED研磨机构对小米进行了研磨实验。由于不能直接测得前述仿真分析得到的碰撞频率和单次碰撞能量，因此本文是通过测量研磨后颗粒的粒径来探究研磨振动频率等参数的影响。研磨后采用Malvern公司制造的Mastersizer 2000激光粒度仪测量颗粒粒度分布。用参数D90(累积体积占总体积的90%的相应粒径)表示研磨效果。f分别为13、 16、 19、 22、 25 Hz，A=25 mm，磨球直径为8 mm，小米初始粒径为1.25~1.43 mm，总质量0.15 g，研磨时间为30 s。
对上述实验所采用的研磨参数通过程序进行研磨仿真。仿真与实验得出的D90如图7所示。D90的仿真结果比实验结果要小，但整体曲线变化趋向基本相同。引起仿真误差的原因可能因为振动频率增大导致碰撞能量迅速增大，使得较细颗粒沉淀和粘附在研磨瓶表面上，阻碍颗粒的破碎进行； 而仿真中，本文假设在任何时刻粉体都被均匀分布在空间中，所以导致仿真结果整体比实验结果小。
|
| 图7 D90随振动频率变化的实验和仿真结果 |
本文采用事件驱动法建立了单个磨球在立式振动球磨机构中的仿真模型。通过该模型分析了振动频率、振动幅度以及研磨罐直径对碰撞频率和碰撞能量两种特征参量的影响。结果表明，在相同的磨腔体积下，随着研磨罐直径从10 mm到14 mm变化，碰撞频率和碰撞能量均减小。碰撞频率、碰撞能量均随振动频率和振动幅度增大而增大，但单次碰撞能量与振动频率成平方关系，而与振动幅度只成线性关系，这表明振动频率比振动幅度对碰撞能量具有更大的影响，在设计振动研磨机构时应优先保证机构具有更高的振动频率。该结论不仅为振动研磨机构设计提供了理论依据，而且为研磨过程的定量研究提供了分析方法。
| [1] | Chen W, Schoenitz M, Ward T, et al. Numerical simulation of mechanical alloying in a shaker mill by discrete element method[J]. KONA Powder Part, 2005, 23:152-162. |
| [2] | Concas A, Lai N, Pisu M, et al. Modelling of comminution processes in Spex Mixer/Mill[J]. Chemical Engineering Science, 2006, 61(11):3746-3760. |
| [3] | Prasad D V N, Theuerkauf J. Improvement in the collision intensity of grinding media in high energy impact mills[J]. Chemical Engineering & Technology, 2010, 33(9):1433-1437. |
| [4] | Ashrafizadeh H, Ashrafizaadeh M. Influence of processing parameters on grinding mechanism in planetary mill by employing discrete element method[J]. Advanced Powder Technology, 2012, 23(6):708-716. |
| [5] | Jayasundara C T, Yang R Y, Yu A B. Effect of the size of media on grinding performance in stirred mills[J]. Minerals Engineering, 2012, 33(3):66-71. |
| [6] | Yamamoto Y, Soda R, Kano J, et al. Dispersion of carbon black by a media mill and correlation of its rate constant with beads impact energy[J]. Powder Technology, 2012, 219(3):105-110. |
| [7] | Pérez-Alonso C, Delgadillo J A. Experimental validation of 2D DEM code by digital image analysis in tumbling mills[J]. Minerals Engineering, 2012, 25(1):20-27. |
| [8] | Gavrilov D, Vinogradov O, Shaw W J D. Simulation of grinding in a shaker ball mill[J]. Powder Technology, 1999, 101(1):63-72. |
| [9] | Pöschel T, Schwager T. Computational Granular Dynamics:Models and Algorithms[M]. Berlin:Springer, 2005. |
| [10] | Gavrilova M L. Empirical studies of optimization techniques in the event-driven simulation of mechanically alloyed materials[J]. The Journal of Supercomputing, 2004, 28(2):165-176. |
| [11] | Reichardt R, Wiechert W. Event driven algorithms applied to a high energy ball mill simulation[J]. Granular Matter, 2007, 9(3):251-266. |
| [12] | Barahona J. Simulation and Optimization of Mechanical Alloying Using the Event-Driven Method[D]. Ottawa, Canada:University of Ottawa, 2011. |
| [13] | Roeller K, Herminghaus S, Hager-Fingerle A. Sinusoidal shaking in event-driven simulations[J]. Computer Physics Communications, 2012, 183(2):251-260. |
| [14] | Richardson D C, Walsh K J, Murdoch N, et al. Numerical simulations of granular dynamics:I. Hard-sphere discrete element method and tests[J]. Icarus, 2011, 212(1):427-437. |
| [15] | Murdoch N, Michel P, Richardson D C, et al. Numerical simulations of granular dynamics:II:Particle dynamics in a shaken granular material[J]. Icarus, 2012, 220(1):296-296. |
| [16] | Schmidt J, Plata M, Tröger S, et al. Production of polymer particles below 5 μm by wet grinding[J]. Powder Technology, 2012, 228:84-90. |