全自动调色机搅拌器叶片流场数值模拟研究

Numerical Simulation of Agitator Blade Flow Field in Fully Automatic Color Mixer

摘要: 以全自动调色机搅拌器为研究对象，采用k-ε湍流模型和多重参考系模型(MRF)，分析不同温度下不同倾斜角度桨叶的流场流动特性，得到搅拌器的应力分布与振动变形，计算不同温度下搅拌器的功率消耗，最后将理论计算和数值模拟结果进行对比。分析结果表明：当搅拌器叶片角度一定时，随着温度的降低，色浆桶内的压力逐渐降低，搅拌器叶片周围的速度范围越来越大；当外界温度一定时，随着叶片倾斜角度的增加，色浆桶内的压力逐渐降低，搅拌器叶片周围的速度范围越来越大，搅拌器的振动可靠性为99.9%，符合工程设计要求，为后续对搅拌器叶片的优化设计提供依据。

Abstract: We take the automatic color mixer as the research object, adopt the k-ε turbulence model and multi-reference system (MRF), analyze the flow characteristics of the blade flow field at different temperatures and different angles, and obtain the deformation of the mixer, the stress distribution and the vibration power consumption of the mixer calculations at different temperatures, and finally the theoretical calculation and numerical simulation results are compared. The analysis results show that when the angle of the agitator blade is fixed, the pressure in the slurry tank gradually decreases with the decrease of temperature, and the speed range around the agitator blade becomes larger and larger. When the external temperature is constant, with the increase of blade tilt angle, the pressure in the color slurry tank gradually decreases, and the speed range around the agitator blade becomes larger and larger. The vibration reliability of the agitator is 99.9%, which meets the engineering design requirements and provides a basis for the subsequent optimization design of the agitator blade.

1. 引言

搅拌设备广泛应用于化工、医药、食品、采矿、涂料、冶金、石油和废水处理等行业中 [1]，其主要作用是促进不同物相之间的传质与传热，因此关于搅拌槽内流体流动状态(包括速度场和压力场的分布与变化)的实验研究和数值模拟一直很受重视 [2]。早期的研究多以实验模拟为主，如选用与真实体系物料性质相似的模拟材料在透明容器中进行搅拌，通过照相技术或直接观察来判断流体运动状态和搅拌效果。随着计算机技术的高速发展，利用计算流体力学对流体运动、传热和相关现象进行分析研究已变得越来越普遍 [3]；全自动色浆长时间放置，会出现凝固或沉淀现象，若是色浆不均匀或者均匀度不够，会严重影响调色精度。在涂料行业 [3]，一般采用搅拌器搅拌色浆，搅拌器设计的好坏将影响色浆的均匀性，进而影响调色精度。使用CFD技术对搅拌槽内流体流动特性和数值模拟研究，可以有效地分析搅拌器是如何实现不同物相之间的传质与传热 [4] [5] [6]，并对后期搅拌器的拓扑优化设计具有重大影响 [7]。运用CFD技术设计搅拌器的同时，通过试验研究搅拌真实流动情况，仍旧是必不可少的过程。党林贵等 [8] 采用计算流体力学仿真软件，研究不同桨叶的搅拌器在搅拌槽内的流动混合特性以及加料位置的变化，结果表明，45˚折叶涡轮桨和平直叶桨的组合类型，可以增强搅拌器内流体的上下湍动，促进混合效果。栾德玉等 [9] 采用流固耦合计算方法，对比分析了错位六弯叶和六弯叶搅拌器的动力学特性，并根据桨叶与流体之间相互耦合运动特性，探讨了宏观流场的结构和搅拌功耗特性，分析了桨叶的变形和等效应力分布。戚振 [10] 基于CFD技术，对搅拌轴单端约束和两端约束的搅拌反应器流体特性进行对比仿真分析，然后根据流固耦合理论对搅拌轴进行有限元分析，观察两种搅拌器的应力应变分布，最后对搅拌反应器以及搅拌器结构振动特性进行分析。综上所述，大多数学者只是通过数值模拟分析搅拌器内流场随着搅拌速度增加所呈现的不同特性，进行流固耦合分析其变形和应力特性分布。很少将外界温度考虑进去，对于粘性液体在不同温度下的搅拌过程不能有效呈现。

本文以全自动调色机搅拌器为研究对象，分析不同温度下不同倾斜角度桨叶的流场流动特性，得到搅拌器的应力分布与振动变形，计算不同温度下搅拌器的功率消耗，最后将理论计算和数值模拟结果进行对比，验证数值模拟结果的合理性，从而得到较为准确的分析结果。

2. 搅拌器的几何模型及相关参数

全自动调色机搅拌器采用浆式搅拌器，其色浆桶内径Φ = 91 mm，搅拌介质(色浆)有效高度为H = 160 mm，搅拌叶片的倾斜角度为17˚、45˚和90˚三种类型，所有类型均采用两层桨叶，共4个桨叶，搅拌桨直径为80 mm，搅拌器叶片宽度为20 mm，底层叶片高度为9 mm，上层叶片高度为69 mm。搅拌器的三维模型图如图1所示。浆式搅拌器带动色浆在搅拌桶内正常工作时，查阅相关资料，得到油性色浆密度为900 kg/m³在常温下的涂四杯粘度为200s，并且温度每降低1℃，其涂四杯粘度增加2~3 s。为了使得设计的搅拌器能在不同温度下对色浆都具良好有搅拌效果，选择在25℃至零下25℃之间，每隔5℃取一个点，10个温度变量下进行设计计算。

Figure 1. 3D model of agitator

图1. 搅拌器三维模型

$P=N_P\rho n{d}^5$ (1)

式中，

P：搅拌功率，W；

N_P：搅拌功率数；

ρ：色浆密度，kg/m³；

n：搅拌器转速，取1 r/s；

d：搅拌桨的直径，m。

永田进治公式是根据搅拌桶内无挡板时，液体流动所形成的“圆柱状回转区”的半径大小与搅拌器叶片所受的流体阻力间的关系推导结果，再通过实验修正所得到的功率特征数计算公式 [11]。虽然该式是针对双桨叶搅拌器而得到的，但是通过对湍流区的实验验证，对于多种搅拌桨叶来说，在桨径相同的条件下，只要桨叶数与桨叶宽度的乘积相等，则所消耗的搅拌功率相等。搅拌作业功率数随流动状态以及搅拌装置的形状和尺寸等条件变化而变化。以往常采用的搅拌功率计算方法有永田进治的关联式计算，对于斜桨无挡板情况下的永田进治计算式如下：

$N_p\text{=}\frac{A}{\mathrm{Re}}+B{\left(\frac{1000+1.2{\mathrm{Re}}^{0.66}}{1000+3.2{\mathrm{Re}}^{0.66}}\right)}^p{\left(\frac{H}{D}\right)}^{\left(0.35\frac{d}{D}\right)}{\sin }^{1.2}\theta$ (2)

$A=14+\left(\frac{b}{D}\right)\left[670{\left(\frac{d}{D}-0.6\right)}^2+185\right]$ (3)

$B={10}^{\left[1.3-4{\left(\frac{b}{D}-0.5\right)}^2-1.14\left(\frac{d}{D}\right)\right]}$ (4)

$p=1.1+4\left(\frac{b}{D}\right)-2.5{\left(\frac{d}{D}-0.5\right)}^2-7\left(\frac{b}{D}\right)$ (5)

$\mathrm{Re}=d^{2}np/\eta$ (6)

式中，

Re：搅拌雷诺数；

θ：搅拌器叶片倾斜角，˚；

b：桨宽，m；

d：搅拌器直径，m；

D：色浆桶直径，m；

η：色浆粘度，Pa∙s；

H：色浆桶高度，m。

3. 搅拌器数值模拟计

以搅拌器三维模型为基础，对其进行单相三维数值模拟分析。选取倾斜角度为17˚、45˚和90˚三种类型搅拌叶片及搅拌轴为流场数值模拟分析对象。

3.1. 搅拌器流场数值模拟值计算

在进行搅拌器数值模拟计算之前，应首先对搅拌流体区域的网格进行离散化，即对于原来的连续区域采用一组有限个离散点进行代替。根据搅拌器流体运动规律，将搅拌槽计算域分为两个部分：一个是转子区域，即包含旋转的搅拌叶片；另一个是定子区域，即包含静止的槽体。整个搅拌器流场数值模拟计算采用多重参考系模型MRF，即转子区域采用旋转参考系，定子区域采用静止参考系，定子区域与转子区域的重合面定义为相互作用面上，假定流动为稳态的，通过相互面来进行流场计算数据的交换 [12] [13] [14] [15]。

在进行网格划分时，网格点之间的邻近关系可分为结构网格、非结构网格和混合网格。网格单元和节点彼此没有固定的规律可循、节点分布完全是任意的是非结构网格的特点；将二者的优势结合起来、同时克服各自的不足是混合网格的优点。本文采用混合网格法对搅拌器进行网格离散化，搅拌器的三维模型与网格划分，如图2所示。

(a) 搅拌器三维模型(b) 搅拌器的网格划分

Figure 2. 3D model and mesh division of agitator

图2. 搅拌器三维模型与网格划分

在对搅拌器流场进行数值模拟计算时，采用基于压力的分离求解器进行计算。其边界条件主要包括：搅拌桶壁面、搅拌轴壁面、搅拌桨叶壁面、动静耦合交界面和动静流动区域设置。其中，搅拌筒体壁面为静态壁面，设置标准静态壁面函数，参数默认设置；搅拌轴壁面和搅拌桨叶壁面为动态壁面，二者的运动分别为：主动运动和被动运动 [16] [17]。因此，搅拌轴壁面设置为绝对运动壁函数，搅拌桨叶运动壁面设置为从动运动壁面函数，并设置相应的搅拌转速。动静耦合交界面设置为interface，用于流场数据的传递与交换；动、静区域均设置为流动区域。搅拌器在不同叶片角度和不同温度下的27种情况下进行流场数值模拟。由于计算结果数据量较大，因此，将结果中的温度为25℃、0℃、−25℃，叶片倾斜角度为17˚、45˚、90˚情况下的压力场和速度场进行对比分析，为使计算结果便于观察，搅拌器压力场分布和速度场分布分别如图3和图4所示。

(a) 25℃下压力场分布 (b) 0℃下压力场分布 (c) −25℃下压力场分布

Figure 3. Pressure field distribution

图3. 压力场分布

由图3搅拌器压力场分布中的9个模型可以看出，最大压力值出现在搅拌器叶片转动方向前端，最小压力值则出现在搅拌器叶片转动方向的后端，且搅拌桶内部的压力整体分布规律是由上部到下部不断减小；在固定水平面，压力值分布则由周边向中心不断减小。当搅拌器叶片角度一致时，随着温度的降低，色浆桶内的压力降低；当外界温度一致时，随着叶片角度的增加，色浆桶内的压力降低。

(a) 25℃下压力场分布 (b) 0℃下压力场分布 (c) −25℃下压力场分布

Figure 4. Velocity field distribution

图4. 速度场分布

由图4搅拌器速度场分布的9个模型中可以看出，搅拌器的速度场分层比较明显，且速度最大值都是出现在搅拌器叶片转动方向前端，且在固定水平面，速度场分布是由周边向中心不断增大。当搅拌器叶片角度一致时，随着温度的降低，搅拌器叶片周围的速度范围越来越大；当外界温度一致时，随着叶片角度的增加，搅拌器叶片周围的速度范围越来越大。速度场的分布情况较为真实反映了搅拌实际状况。

3.2. 搅拌器流固耦合分析

搅拌器的材料选择ABS塑料，搅拌器的叶片角度选择45˚，在−25℃的温度下进行分析。利用有限元分析软件的接口将得到的流固耦合交界面的压力数据传递到搅拌器结构的表面，进行耦合场计算。在搅拌器底部添加固定约束，搅拌器转速为50 r/min。搅拌器网格划分模型，如图5所示，整个搅拌器的单元数为12,369，节点数为23,375。

在流体压力与离心力耦合作用下，进行搅拌器的变形和应力分析，得到搅拌器位移变形和应力分布云图，由图6和图7可知，搅拌器的变形最大位置在搅拌器叶尖部位，最大变形量为0.030073 mm，属于合理变形；应力最大位置在叶片轴根部位，最大应力为0.086354 MPa，符合强度要求。搅拌器下层叶片变形量比上层变形量大，下层叶片的应力值比上层叶片的应力值小，主要是因为底部是封闭状态，流体运动速度快，故变形较大应力较小。变形和应力的分布情况较为真实反映了搅拌实际状况。

Figure 5. Mixer meshing model

图5. 搅拌器网格划分模型

将静力学分析得到的数据作为预应力传递到谐响应分析中，对搅拌器进行谐响应分析，得到搅拌器的振动变形与应力分布。由图8和图9可知，搅拌器的最大振动变形位于搅拌器底部的叶片处，最大振动应力位于搅拌器底部搅拌轴的根部，比较符合实际搅拌情况。下方由于节流孔的存在，出现漩涡。

Figure 6. Agitator displacement deformation

图6. 搅拌桨振动等效应力

Figure 7. Equivalent stress of agitator

图7. 搅拌桨等效应力

Figure 8. Agitator vibration displacement deformation

图8. 搅拌桨位移变形

Figure 9. Equivalent stress of agitator vibration

图9. 搅拌桨振动位移变形

3.3. 搅拌器可靠性分析

搅拌器设计完成之后，为了验证搅拌器结构设计参数(如：搅拌器结构、搅拌器材料、搅拌器转速等)的准确性，对其进行可靠性分析。目前对可靠性分析方法 [18] [19] [20] 很多，如蒙特卡罗法、一次二阶矩阵法、故障树法、响应面法等，其中，响应面法在提高计算效率的同时能满足计算精度，故本节选用响应面法对搅拌器振动进行可靠性分析。为了探究在不确定因素影响下搅拌器振动可靠性，选择搅拌器材料密度、弹性模量，转速作为随机输入变量，各变量间相互独立且服从正态分布，随机变量统计特征，如表1所示。搅拌器最大振动变形及最大振动应力作为输出响应，查阅相关资料可以得到，搅拌器许用振动变形为8.0 × 10⁻⁶ mm，许用振动应力为7.2 × 10⁻⁵ MPa。

利用拉丁超立方抽样对随机变量进行100组抽样，并代入有限元分析模型中分析得到各组样本对应的输出响应，利用抽取的样本建立响应面方程，如式(7)所示，用该方程代替搅拌机振动有限元分析过程，与蒙特卡罗法相结合进行大批量抽样，对抽样结果进行统计分析，得到搅拌器振动变形与应力的分布直方图，如图10所示。由图10可看出，在不确定因素影响下，搅拌器的振动变形与应力大体服从正态分布。将搅拌器振动变形与振动应力的抽样结果与搅拌器许用振动变形与许用振动应力比较，得出在流固耦合影响下，搅拌器的振动可靠性为99.9%，符合工程设计要求。

Table 1. Statistical characteristics of random variables

表1. 随机变量统计特征

$\{\begin{array}{l}{y}_1=-0.016+0.003x_1+8.502x_2+0.014x_3-0.002x_1^2-2.096x_2^2-0.008x_3^2\\ y_2=-0.012+0.025x_1+6.639x_2+0.115x_3-0.013x_1^2-1.636x_2^2-0.068x_3^2\end{array}$ (7)

式中：y₁为搅拌器振动变形；y₂为搅拌器振动应力；x₁为搅拌器材料密度；x₂为搅拌器材料弹性模量；x₃为搅拌器转速。

(a) 振动变形直方图(b) 振动应力直方图

Figure 10. Histogram of vibration deformation and stress distribution of agitator

图10. 搅拌器振动变形与应力的分布直方图

4. 理论计算和数值模拟

不同温度下搅拌器消耗功率理论计算的结果和流场数值模拟结果进行对比分析，验证流场数值模拟分析的合理性。

通过前期对搅拌器理论计算，可以得到其理论消耗功率，如表2所示。在对搅拌器进行流体模拟计算时，可以得到搅拌器搅拌油漆的力矩数值，并且与相应搅拌转速相乘可计算出流场数值模拟中搅拌器实际消耗功率，如表3所示。

将表2、表3中理论计算数据与数值模拟计算数据进行比较，可知理论计算数据与数值模拟计算数据较为接近，吻合度良好，为便于更直观的表示，绘制两者实际消耗功率曲线如图11所示。

由图11可知，理论计算曲线趋势与数值模拟曲线趋势一致，数值模拟值略高于理论计算值，其原因为理论计算认为搅拌桨为无摩擦旋转运动，没有将摩擦阻力考虑在其中，所以造成理论计算值低于数值模拟值。由表2可得到，在室温25℃下，搅拌器理论消耗功率为0.0065 W；由表3得到，在室温25℃下，搅拌器数值模拟消耗功率为0.0075 W。则搅拌器消耗功率理论计算值、数值模拟分析值偏差在15%左右，符合要求。既验证在25℃温度下，搅拌油漆消耗功率的理论计算结果和流场数值模拟结果的正确性，又进一步验证了流场数值模拟计算方法的合理性。

Table 2. Agitator theoretical power consumption

表2. 搅拌器理论消耗功率W

Table 3. Agitator numerical simulation power consumption

表3. 搅拌器数值模拟消耗功率W

(a) 叶片角度17˚ (b) 叶片角度45˚(c) 叶片角度90˚

Figure 11. Mixer theory and numerical simulation comparison chart

图11. 搅拌器理论与数值模拟对比图

5. 结语

1) 压力场分析结果：在固定水平面，压力值分布则由周边向中心不断减小。当搅拌器叶片角度一致时，随着温度的降低，色浆桶内的压力逐渐降低，搅拌器叶片周围的速度范围越来越大；

2) 速度场分析结果：在固定水平面上，速度场分布是由周边向中心不断增大。当外界温度一致时，随着叶片倾斜角度的增加，色浆桶内的压力逐渐降低，搅拌器叶片周围的速度范围越来越大；

3) 搅拌器的变形最大位置在搅拌器叶尖部位，搅拌器下层叶片变形量比上层变形量大，下层叶片的应力值比上层叶片的应力值小；

4) 搅拌器的振动变形与应力大体服从正态分布，搅拌器的振动可靠性为99.9%，符合工程设计要求。在25℃室温下，其实际消耗功率和理论计算值与数值模拟值误差范围在20%以内，符合要求。