搅拌功率计算程序的开发

王智勇, 秦 泰, 刘 刈, 罗云蓉, 唐克伦

(1.农业部沼气科学研究所, 成都 610041;2.成都市城市道路桥梁管理处, 成都 610017;3.四川理工学院机械工程学院, 四川 自贡 643000)

搅拌功率计算程序的开发

王智勇1, 秦 泰2, 刘 刈1, 罗云蓉3, 唐克伦3

(1.农业部沼气科学研究所, 成都 610041;2.成都市城市道路桥梁管理处, 成都 610017;3.四川理工学院机械工程学院, 四川 自贡 643000)

搅拌设备的设计除应保证工艺要求外,以较小功率消耗达到搅拌目的亦是设计目标之一。搅拌器功率计算一方面解决了需要向搅拌介质提供多少功率问题，以便有效地选择电机；另一方面为搅拌器强度计算提供依据，以保证桨叶、搅拌轴的强度。因而，搅拌器功率计算是搅拌器设计中必不可少的步骤。然而搅拌器功率计算却非常繁琐。将Rushton图算法中算图曲线数据进行人工采集，利用仿射变换对采集的数据进行了修正，最后利用B样条曲线对Rushton算图参数化插值，并利用Rushton图算法与永田进治公式法的互补特性，完成了对常规叶片形状的搅拌器功率计算程序化，以达到减小搅拌器功率计算的工作量的目的。

搅拌器;功率计算;Rushton算图;永田进治公式;程序设计

引 言

使两种或多种不同的物质在彼此之间互相分散，从而达到均匀混合的单元操作称为物料的搅拌[1-2]。搅拌是化工生产过程常见的单元操作装置中的重要组成部分，通过搅拌可以使互溶物流均匀溶解，使不互溶的物料很好地分散或悬浮，强化了传热或传质过程。正因为搅拌操作具有上述功能，搅拌设备在工业诸多领域中有着广泛的用途，特别是在制药、精细化工、染料、油漆等产品的生产中应用更为广泛。

在满足所需工艺参数的前提下，搅拌设计的目标是利用最小的功率消耗达到搅拌的目的。对于不同的搅拌目的和搅拌过程，所需要的搅拌功率是不同的。在固-液悬浮操作中，有时需要通过搅拌来防止固体颗粒发生沉降，这时只需搅拌器输入较小的能量，使全体固体颗粒在搅拌槽底部浮游起来即可，并不要求固体颗粒均匀地分布于整个槽体内；而对于均匀化搅拌，为了达到搅拌均匀程度的要求，根据搅拌的工艺参数以及搅拌的转速等要求，需要消耗的功率将增大[2-3]。

不同的搅拌目的对于搅拌功率的要求不尽相同。计算搅拌功率的目的有二：一是为了解决一定型式的搅拌器能向被搅拌介质提供多大功率的问题，以满足搅拌过程的要求，并选配合适的电机；二是为搅拌器的强度计算提供依据，以保证桨叶、搅拌轴的强度[1]。正确地计算搅拌作业所需的功率对于节约能量和提高搅拌操作效果非常重要。在传统的搅拌设备设计过程中，设计人员必须采用手工计算的方式测算搅拌功率，由于其涉及到的计算多，需要花费较长的时间[2-3]，进而影响了设计过程。随着计算机技术高速发展，搅拌器内流体仿真逐渐成为可能。为了对搅拌器流场有深入了解，同时也能够顺便计算出搅拌器功率，越来越多用户选择了基于CFD的流场仿真分析[6-11]。基于CFD流场分析能够实现各种形式、各种工况(包括多相流)[12-13]、各种安装方式的搅拌器流场分析以及其轴功率计算。不过，基于CFD流场分析计算功率对设计者要求较高，需要掌握流体力学相关专业知识以及专业的流体仿真软件；同时，进行CFD流场分析[14-17]，其建模过程、流体动力学分析过程也是非常耗时的，流体动力学分析过程动辄需要连续计算几十乃至几百上千小时。而且，几十万元人民币的专业CFD软件，加上高档计算机配置对于一般公司也是不可小觑的消费。

为了简化设计过程，为设计人员提供参考，本文着重探讨将常用的搅拌功率计算程序化的过程。

1 影响搅拌功率的因素

搅拌器的功率与搅拌槽内流体的流动状态有关。影响搅拌功率的因素很多，其中包括结构的影响，如搅拌桨的直径、宽度、倾斜角度、转速；搅拌槽参数的影响，如槽体的直径、深度、有无挡板等条件；被搅拌物料特征的影响，如密度、粘度等[1-3]。虽然影响搅拌功率的因素很多，这些因素归纳起来可以称之为桨、槽的几何参数，桨的操作参数以及影响功率的物性参数。

2 计算原理

2.1 量纲分析

在分析搅拌功率的过程中，为简化分析过程，可以假定桨、槽的几何参数均与搅拌器的直径有一定的比例关系，并将这些比值称为形状因子。对于特定尺寸的系统，形状因子一般为定值，故桨、槽的几何参数仅考虑搅拌器的直径d。桨的操作参数主要指搅拌器的转速n，物性参数主要包括搅拌流体的密度ρ和黏度μ。但搅拌发生打旋现象时，重力加速度g也将影响到搅拌功率。因此搅拌功率与各变量之间的关系可以表示为[1]：

P=f(n,d,ρ,μ,g)

(1)

也可以写成指数形式：

P=Kna1da2ρa3μa4ga5

(2)

在式(1)与式(2)中，P为搅拌功率，单位：W；K为系数，与系统的几何形状有关；n为搅拌转速，单位：r/s；d为搅拌器的直径，单位：m；ρ为流体的密度，单位：kg/m3；μ为流体的黏度，单位：Pa.s；g为重力加速度，单位：m/s2；a1～a5为待定常数。对式(2)进行量纲分析可得：

(3)

即

(4)

令Np=P/ρn3d5，称为功率特征数；令Re=ρnd2/μ，称为搅拌雷诺数，表示流体惯性力与粘滞力之比，用于衡量流体的流动状态；令Fr=n2d/g，称为弗劳德数，表示流体惯性与重力之比，用以衡量惯性的影响；令x=-a4，y=-a5。则式(4)可以写为：

Np=KRexFry

(5)

若再令φ=Np/Fry，称为功率因素，则有：

φ=KRex

(6)

从量纲分析法得到搅拌功率特征数关联式后，可对一定形状的搅拌器进行一系列的实验，找出各流通范围内具体的经验公式或关系算图，则可以解决搅拌功率的计算问题。

2.2 搅拌功率的计算

结合量纲分析法，研究人员通过研究均相系统的搅拌得到了其搅拌功率的计算表达式，并以此为基础达到了研究非均相系统搅拌功率的目的。

均相物系搅拌功率计算有很多种方法，常用的是Rushton图算法和永田进治公式法。

2.2.1 Rushton图算法

Rushton算图适合于推进式、涡轮式和桨式搅拌器搅拌功率的计算[1-3]。Rushton等人对多种型式的搅拌器在液体黏度为1×10-3～40 Pa.s、Re<106时进行了大量实验，测定了各种条件下的搅拌功率，并整理得出功率因素φ和雷诺数Re的关系图，如图1所示。图1中纵坐标为功率因素φ，横坐标为Re。由图1中的曲线可以看出：搅拌槽中的流体的流动可根据Re的大小大致分为三个区域，即层流区、过渡区和湍流区。

由图1可知，在低雷诺数(Re≤10)的层流区，流体不会打旋，由于反映重力影响的Fr可以忽略，φ=Np。所以在层流区，不同搅拌器的功率因素φ为斜率为-1的直线；当10≤Re≤10 000时为过渡区，功率因素为一下凹曲线；当Re>10 000时，流体进入充分湍流区，为了消除打旋现象，一般采用全挡板条件，故重力的影响可以忽略不计，所以φ=Np，此时功率曲线呈现水平直线状态，即功率因素φ与Re无关，为一恒定值。

搅拌功率特征数Np是搅拌设备最基本的特征参数之一。搅拌功率按下式计算[1]：

P=Npρn3d5

(7)

由于Rushton算图对搅拌器的类型和搅拌器的几何比例关系有严格的限制，使用Rushton算图查对应的功率因素φ时，如果比例关系不同，功率因素φ也不同。

2.2.2 永田进治公式法

日本永田进治等根据在无挡板直立圆槽中搅拌时“圆柱状回转区”半径的大小及桨叶所受的流体阻力进行理论推导，并结合实验结果确定了一些系数而得出双叶搅拌器功率的计算公式。永田进治公式可近似用于桨式、多叶开启涡轮、圆盘涡轮等常用桨型无挡板湍流区搅拌功率的计算[1]。永田进治公式为：

图1 Rushton：φ-Re关系算图

(8)

式中：

b—搅拌桨叶的宽度,m；

H-液层的深度，m；

D-搅拌槽的内径，m；

d—搅拌器的桨径,m；

θ—桨叶的折页角，对于平桨θ=90度。

2.2.3 图算法与公式法的关系

在常用的推进式、涡轮式、开启式涡轮、闭式涡轮和桨式搅拌器等常用的搅拌功率计算过程中，使用Rushton算图可以计算推进式、涡轮式和桨式搅拌器的搅拌功率，其中计算三叶推进式和六片平直叶圆盘搅拌器搅拌功率时，搅拌器可以有挡板也可以无挡板；在计算六片弯叶圆盘涡轮、六片箭叶圆盘涡轮等搅拌器的功率时，搅拌器必须有挡板；而永田进治公式可以近似计算桨式、多叶开启式涡轮、圆盘涡轮等常用桨型无挡板区搅拌功率的计算。因而，程序设计采用Rushton图算法和永田进治公式法计算搅拌功率可以比较好地形成计算互补，见表1。

表1 图算法与公式法的适用范围

注：s-桨叶螺距，d-搅拌器直径，Y-有挡板，N-无挡板。

3 程序的开发

Microsoft的Visual Basic 6.0是基于Windows平台的应用程序开发工具。Visual Basic提供了完善的可视化编程环境，程序的运行基于事件驱动机制，大大减低了程序的编写难度。程序设计人员可以快速设计界面、编写程序、调试程序直至把应用程序编译输出成文件后缀名为“.EXE”的可执行文件[4]。使用Visual Basic可以有效地提高应用程的运行效率和可靠性。

程序开始执行后，如果用户选择“读入最近一次的计算数据”，程序会从保存数据的文件中读取最近一次计算使用的数据并以此为初始值填入各输入窗口。输入数据结束后，程序会检测输入数据的完整性，并依次计算雷诺数Re和挡板系数；接下来，程序会根据用户输入的搅拌器的类型和挡板情况选择使用Rushton图算法或是永田进治公式法来计算搅拌功率。程序的流程如图2所示。

图2 程序流程

4 程序特点及功能实现

(1) 通过Rushton算图获得功率因素φ

Rushton算图是在大量实验的基础上整理得到的(图1)，图1中的纵坐标为功率因素φ，横坐标为雷诺数Re，共有8种桨型的搅拌器在有挡板和无挡板条件下的12条关系曲线。使用Rushton算图计算搅拌功率，首先要根据相关的数据计算得到雷诺数Re，根据得到的Re数在Rushton算图中查得功率因素φ进而可以计算得到搅拌功率。在整个计算过程中，查表非常关键。如果查表得到的数据准确，就可以得到相对精确的计算结果，如果查表精度不高，则计算误差会增大。在图1中，12条曲线的Re数从小于10跨越到大于104,不可能通过计算所有的Re数值来建立检索列表供查询时使用。因此如何选择计算方法从而能快速准确地得到功率因素是重点研究内容。

分析图1中的12条曲线，所有的曲线均表现光滑，且不存在断点，因此如果将12条曲线采用数据逼近的方式逼近之后，再通过插值计算得到所需结果，可以成为解决问题的方法之一。基于此思想，本文获取功率因素φ的流程如下：

步骤1 获取Rushton算图中的数据点。在Matlab中读入并显示Rushton曲线，利用Matlab中“data cursor”[18]工具在每条曲线上均匀采样获取每条曲线的数据点坐标(像素坐标)。在Matlab中，图像可以进行缩放，用“data cursor”工具在放大的图像上取点时，可以很方便地点选到曲线上的像素点的准确坐标。使用该方法获得的关键点的数据比人工读取时精度更高，数据更准确。

步骤2 图像的仿射变换。考虑到所采集的图形可能带有部分倾斜和平移，在Matlab中利用“data cursor”工具所获取的数据点不能正确反应极坐标的数据值；同时，Matlab采集的坐标为图形的像素坐标，需要转换到相应的Rushton算图坐标系下的坐标。因而，需对图形作相应的仿射变换，求出两种坐标系下数据点的关系。其变换矩阵[19]关系为：

(9)

用“data cursor”工具获得Rushton算图坐标系的矩形四个角点像素坐标，其相应的数据点坐标也已知。这样通过式(9)建立8个方程，用最小二乘法即可求出旋转和缩放矩阵R与平移向量T中的各元素。由此，步骤1中获得的数据点像素坐标值即可通过式(9)转换到相应的Rushton算图坐标系下数据点坐标值。

步骤3 B样条曲线逼近。B样条是一个称为节矢量的非递减的参数u的序列U：u0≤u1≤…≤ui+k+1所决定的k次分段多项式，也即是k次多项式样条[20]。其中三次B样条曲线C2阶连续，用于数据逼近可以获得连续光滑的曲线。针对前面采集的每条曲线上约230个数据点，采用三次B样条曲线逼近。通过逼近尝试，当曲线段为10时，即已较好地逼近了原始曲线。综合考虑，本文采用了12条曲线段实施曲线的逼近。对每一条逼近曲线，保存其相应的控制点。

步骤4 在Visual Basic 6.0中，用步骤3保存的控制点得到B样条曲线，通过二分法计算雷诺数Re所对应的功率因素φ。

在上述的过程中，步骤1～步骤3需要借助Matlab实现。在搅拌功率计算的过程中，只有步骤4是在VB中实现的。通过Matlab处理之后，复杂的Rushton算图通过B样条控制点反算得到B样条曲线即可再现算图中的曲线，比较图1和图3，两图曲线基本一致。同时分析图3的12条曲线，曲线表现光滑，没有出现尖点和断点，因此可以说明文章采用的方法是可行的，得到的功率因素相对也是准确可靠的。

图3 三次B样条曲线逼近的Rushton算图

(2) 对比研究

由于计算搅拌功率的目的之一就是为了解决一定型式的搅拌器能向被搅拌介质提供多大功率的问题，以满足搅拌过程的要求，并选配合适的电机。当某一次的功率计算完成之后，用户通常希望只修改少数几个参数就可以重新计算得到搅拌功率，以达到对比的目的。在程序中，如果用户在开始新计算时，通过勾选“读入最近一次的输入数据”选项，程序可以自动载入最近一次的相关计算数据，用户在此基础上修改少量的数据即可获得新的搅拌功率，实现对比研究的目的，从而为设计者确定桨、槽的几何参数和选择合适的电机作参考。

记录并保存最近一次的相关数据是实现对比研究的基础。记录和保存用户输入的数据，可以在三个时期完成：①在用户输入数据的同时进行保存；②计算搅拌功率之前进行保存；③退出程序前保存数据。如果采取在用户输入数据的过程中保存数据，降低了的程序的可操作性，减缓了程序的执行，同时也不利于程序的编写。如果在退出程序前保存数据，有可能保存的数据并非计算中所使用的数据，造成对比研究的失败。选择在计算搅拌功率前保存数据，则可以保证记录的数据与当前计算所使用的数据一致，同时程序的编写也比较简单。综合考虑上述三种方案，在计算搅拌功率前保存数据是合适的，也是合理的。软件界面如图4所示。

图4 搅拌器功率计算软件界面

5 结束语

利用Visual Basic 6.0编程，将计算搅拌功率的计算过程程序化是一个有益的尝试。在设计搅拌装置的过程中，通过文章设计的计算程序，可以方便而迅速地依据搅拌槽、桨的几何参数，搅拌液的物理特性等计算得到搅拌器转轴的功率。在设计搅拌设备时，通过理论计算与设计者的实际经验相结合，可以为搅拌设备的设计带来便利。

实际的搅拌器叶片种类多种多样，并且处于不断开发之中[5]。文章主要采用Rushton图算法和永田进治公式法，是目前应用较多的两种算法，包含了目前广泛应用的几种搅拌器叶片形式，算法上两者刚好实现互补性。该程序能够解决目前搅拌器设计中经常遇到的搅拌器叶片的搅拌器功率计算问题，为搅拌设备的设计提供极大的便利。

[1] 柴诚敬.化工流体流动与传热[M].北京:化学工业出版社,2007.

[2] 王凯,冯连芳.混合设备设计[M].北京:机械工业出版社,2000.

[3] 陈志平,章序文,林兴华.搅拌与混合设备设计选用手册[M].北京:化学工业出版社,2004.

[4] 江建卫,冯忠绪,毛兴中.搅拌功率的计算[J].西安公路交通大学学报,2000,20(3):98-100.

[5] 陈志希,谢明辉,周国忠,等.24种搅拌器的功率曲线[J].化学工程,2010,38(3):38-43.

[6] 梁敬福,黄振峰,李欣欣.大型沼气厌氧池侧进式搅拌流场的数值模拟分析[J].化学反应工程与工艺,2017,33(1):82-89.

[7] 方德明,陈涛,杨象岳,等.基于CFD流场分析的反应釜搅拌器结构改进[J].轻工机械,2014,32(2):95-98.

[8] 肖伟,王凡宇.基于Fluent的双卧轴搅拌器数值模拟研究[J].机械研究与应用,2016(4):39-41.

[9] 丁健华,马腾,陈涛,等.搅拌釜内流场三维数值模拟及功率预测[J].化工装备技术,2015(1):19-22.

[10] 党林贵,郭淑雪,王定标,等.不同组合桨搅拌器搅拌特性的数值研究[J].郑州大学学报:工学版,2013,34(3):59-62.

[11] 谭妍玮,张众华,刘秋生.基于Fluent的考虑旋转车轮影响的汽车外流场分析[J].四川理工学院学报:自然科学版,2015,28(4):17-21.

[12] 唐克伦,张应迁,梁智权.多相流搅拌器流场数值模拟软件[J].计算机辅助工程,2011,20(1):37-41.

[13] 白晓莉,唐克伦,李辉,等.基于Fluent的搅拌器三维流场数值模拟及其实验研究[J].机械工程师,2015(4):17-21.

[14] 徐建航,杨翊仁.导向管底部泄流孔水力损失计算分析[J].四川理工学院学报:自然科学版,2016,29(6):44-47.

[15] 李照广,张全贞,吕建.钢结构焊接的Ansys数值模拟[J].四川理工学院学报:自然科学版,2016,29(5):38-42.

[16] 徐伟,杨翊仁,文华斌.高速列车通过隧道时动态数值模拟[J].四川理工学院学报:自然科学版,2015,28(5):16-20.

[17] 付磊,付丽娅,唐克伦,等.基于FLUENT的管壳式换热器壳程流场数值模拟研究[J].四川理工学院学报:自然科学版,2012,25(3):17-21.

[18] 陈杰.MATLAB宝典[M].3版.北京:电子工业出版社,2011.

[19] 潘云鹤.计算机图形学:原理、方法及应用[M].北京:高等教育出版社,2002.

[20] 朱心雄.自由曲线曲面造型技术[M].北京:科学出版社,2000.

Programming of the Shaft Power Computing of the Stirrer

WANGZhiyong2,QINTai1,LIUYi1,LUOYunrong3,TANGKelun3

(1.Biogas Research Institute of the Ministry of Agriculture,Chengdu 610041,China;2.Municipal Facility Management Office of Chengdu, Chengdu 646000,China;3.School of Mechanical and Engineering, Sichuan University of Science & Engineering,Zigong 643000,China)

The design of stirring equipment should not only guarantee the technological requirements, but also achieve a less power consumption to meet the stirring goals. The agitator power calculation has solved the problem of how much power need to provide to the mixing medium, so as to select the motor effectively; on the other hand, the basis is provided for the stirrer strength calculation, to ensure strength of the blade and the stirring shaft. Therefore, stirring power calculation is essential in the design process. However the agitator stirring power calculation is very complicated. The nomogram curve data of Rushton graph algorithm are artificially collected, then the collected data are corrected by affine transformation, and finally are parameterized and interpolated by the B-spline curve. Due to the complementary characteristics of Rushton graph algorithm and Nagata Shinji formula method, the mixing power on conventional vane shape is programmed to lessen the design workload.

stirrer; power computing; Rushton graph; Nagata Shinji formula; programming

2017-06-22

国家重点研发计划课题(2016VFD0501403);材料腐蚀与防护四川省重点实验室基金(2016CL17)

刘 刈(1983-),男,辽宁鞍山人，助理研究员,博士，主要从事畜粪污沼气化处理利用的研究,(E-mail)Liuyi@caas.co

1673-1549(2017)04-0017-06

10.11863/j.suse.2017.04.04

TP311.1

A