2. 绿色制药技术与装备教育部重点实验室, 浙江 杭州 310023;
3. 浙江工业大学 机械工程学院, 浙江 杭州 310023
2. Key Laboratory for Green Pharmaceutical Technologies and Related Equipment of Ministry of Education, Zhejiang University of Technology, Hangzhou 310023, China;
3. College of Mechanical Engineering, Zhejiang University of Technology, Hangzhou 310023, China
在搅拌反应器的工程设计和运行调控过程中，反应器搅拌功率的准确预测是一个较为重要的因素[1]。采用传统经验公式预测搅拌功率，由于经验公式的适用范围和对象有限，其准确度通常不高，且无法完全满足工程实际的需要。随着计算流体力学(CFD)的发展，工程设计中采用CFD软件进行搅拌功率预测，对于纯液体系而言，预测效果总体表现较好[2-3]。然而，对于化工、冶金、食品和制药等过程工业中需要借助搅拌来强化的固-液两相体系，特别是固含率大于20%的高固稠密体系[4-5]，运用CFD方法预测搅拌功率的相关研究大都处于定性层面[6-7]，搅拌功率的CFD预测值与实际值之间的偏差情况并不明确。因此，在实际应用中，通过CFD方法来预测固-液体系的搅拌功率缺乏一定的可信度。
Gu等[8]采用Euler-Euler方法和标准k-ε湍流模型对水-玻璃珠体系进行了模拟，得到了3种不同桨叶功率随转速的变化曲线。李晨辰[9]采用RNG k-ε模型与Euler-Euler方法研究了固含率为17%的水-石英砂体系的搅拌功率，并对比了不同转速下的实验结果与模拟结果，发现平均偏差为9. 8%。总体来看，对于固-液两相体系搅拌功率的研究，Euler-Euler方法的应用较为广泛。Euler-Euler方法中常用于固-液两相流模拟的Euler模型与Mixture模型，是将固相当量黏度与液相黏度(或是两者黏度与剪切速率的乘积)按照各自的体积分数进行加权，得到固-液两相体系的黏度和黏性应力。Konijn等[10]对悬浮液黏度的实验研究结果表明，悬浮液的黏度随着固含率的增大有着较明显的增大。而按软件中基于颗粒动力学得到的固相当量黏度，相对于液相黏度而言非常小，导致固含率的增大不会引起悬浮液体系黏度的显著增大。因此，Euler模型与Mixture模型直接用于固-液搅拌功率的预测并不合适，有必要对其中的黏度计算模型进行修正。
为此，本研究以甘油-聚甲基丙烯酸甲酯(PMMA)悬浮液体系为研究对象，基于实验结果提出了考虑固含率和剪切速率的悬浮液体系黏度模型，并结合商用Fluent软件对固含率为0. 1~0. 4工况下的固-液两相流搅拌釜搅拌功率进行了数值模拟，探究了该模型在不同固含率及不同搅拌转速的适用性。此外，模拟结果亦与实验结果进行对比验证。
2 实验及数值模拟 2.1 实验装置实验装置如图 1所示，搅拌釜为带挡板的有机玻璃椭圆底圆柱形搅拌釜，搅拌釜内径为0. 150 m，高H1为0. 220 m，液面高度H2为0. 190 m(此液面高度下总体积为3. 27 L)，挡板宽度为0. 013 m，厚度为0. 002 m。搅拌桨的直径为0. 075 m，桨叶距离釜底高度H3为0. 122 m，搅拌桨搅拌方向在俯视下为顺时针方向。装置中的扭矩传感器型号为DYN-200，量程为±0. 1 N ·m，精度为0. 1%。
|
图 1 搅拌实验装置 Fig.1 Schematic diagram of the mixing setup 1-motor 2-frequency modulator 3-data acquisition 4-computer 5-thermometer 6-thermocouple 7-blade 8-stirring tank 9-baffle 10-torque sensor |
实验原料采用的液相介质是分析纯甘油，其密度为1 238 kg ·m-3，黏度为0. 567 Pa ·s (30. 5 ℃，本研究中黏度指动力黏度)，固相介质为PMMA颗粒，颗粒直径为0. 150 mm，密度为948. 2 kg ·m-3。
2.2.2 搅拌功率测量方法由于实验研究的对象涉及高固含率(文中固含率皆为体积固含率)固-液两相体系，所需加入的PMMA固体颗粒所占体积较大，不论桨叶处在较高位置还是较低位置，搅拌后都无法使颗粒完全进入液相，影响后续的实验研究。为解决这一问题，选择在加入固体颗粒后先用玻璃棒进行预搅拌，使得固体颗粒全部润湿。静置一段时间后，由于PMMA固体颗粒密度小于甘油密度，固体颗粒会聚集在搅拌介质最上方形成一定厚度的固-液混合层，而下方则为不含固体颗粒的清液层。经测量，当总体固含率在0. 1~0. 4时，混合层固含率均为0. 55，以此作为后续实验以及模拟的初始条件。此外，考虑甘油的黏度对温度比较敏感，而搅拌过程会使釜内温度升高，因此，在距釜底40 mm与150 mm的侧壁处，有两根连接着测温仪的热电偶用于监测釜内介质的温度。通过监测发现整个实验过程中搅拌介质的温度变化在30~31 ℃，因此可忽略温度变化造成的影响。
搅拌功率测量实验涉及7个不同转速和4个不同固含率，总共28种工况，具体的转速和固含率如表 1和表 2所示。
|
表 1 搅拌实验转速汇总 Table 1 Summary of solid contents in the mixing experiments |
|
|
表 2 搅拌实验固含率汇总 Table 2 Summary of solid contents in the mixing experiments |
搅拌功率采用扭矩测量法进行测量，实际搅拌扭矩为负载扭矩MF减去空载扭矩MK，再根据式(1)计算搅拌功率P。
| $ P=2 \pi N\left(M_{\mathrm{F}}-M_{\mathrm{K}}\right) $ | (1) |
式中：P为搅拌功率，W；N为转速，r ·min-1；MF为负载扭矩，N ·M；MK为空载扭矩，N ·M。
2.2.3 黏度测量方法甘油与PMMA固体颗粒按不同体积比混合得到不同固含率的悬浮液体系，充分搅拌使其中固体颗粒尽可能分布均匀。使用哈克旋转流变仪(HAAKE Rheometer，赛默飞)测量不同固含率的悬浮液体系在不同剪切速率下的黏度，剪切速率变化范围在1~300 s-1。由于搅拌功率测量实验中的介质温度在30~31℃变化，黏度测量实验过程中悬浮液的温度控制为30. 5 ℃。
旋转流变仪通过已知的转子剪切速率γ以及所测得的剪切应力τ按式(2)计算得到动力黏度值。
| $ \mu=\frac{\tau}{\gamma} $ | (2) |
式中：μ为动力黏度，Pa ·s；τ为剪切应力，Pa；γ为转子剪切速率，s-1。
2.2.4 取样点选取及固含率测量方法选择距釜底65、100、135、170 mm高度处的4个排液口为取样点，从排液口处取得样品后称得质量为m，通过真空泵、漏斗和滤纸过滤得到样品中的固体颗粒，再通过鼓风干燥箱干燥后测量得到固体颗粒的质量ms，最后结合固相密度ρs和液相密度ρl按式(3)计算得到固含率φ。
| $ \varphi=\frac{\frac{m_{\mathrm{s}}}{\rho_{\mathrm{s}}}}{\frac{m_{\mathrm{s}}}{\rho_{\mathrm{s}}}+\frac{m-m_{\mathrm{s}}}{\rho_{\mathrm{l}}}} $ | (3) |
式中：φ为固含率；ms为干燥后固体颗粒的质量，kg；ρs为固相密度，kg ·m-3；m为样品质量，kg；ρl为液相密度，kg ·m-3。
2.3 模拟计算方法本研究分别使用Euler模型、Mixture模型以及对固-液两相体系黏度计算模型进行修正后的Mixture模型对不同固含率下的固-液搅拌功率进行了模拟。
2.3.1 固-液两相体系黏度计算公式Euler模型与Mixture模型这两种模型对于固-液两相体系中黏性应力的计算式如式(4)~(5)。
| $ \begin{aligned} &\text { Euler 模型 } &\tau_{\mathrm{m}}=\varphi_{\mathrm{s}} \mu_{\mathrm{s}} \gamma_{\mathrm{s}}+\varphi_1 \mu_{\mathrm{l}} \gamma_1 \end{aligned} $ | (4) |
| $ \text { Mixture 模型 } \quad \tau_{\mathrm{m}}=\left(\varphi_{\mathrm{s}} \mu_{\mathrm{s}}+\varphi_1 \mu_{\mathrm{l}}\right) \gamma_{\mathrm{m}} $ | (5) |
式中：τm为混合液的黏性应力，Pa；φs为体积固含率；φl为液相体积分数；μ为动力黏度，Pa ·s；γ为剪切速率，s-1；下标s表示固相；下标l表示液相；下标m表示混合液。
2.3.2 对Mixture模型的黏度修正方法Fluent软件中对于固相当量黏度μs的计算如式(6)。
| $ \mu_{\mathrm{s}}=\mu_{\mathrm{s}, \mathrm{col}}+\mu_{\mathrm{s}, \mathrm{skin}}+\mu_{\mathrm{s}, \mathrm{fr}} $ | (6) |
式中：μs为当量黏度，Pa ·s；μs, col为固相碰撞黏度，Pa ·s；μs, skin为固相动力黏度，Pa ·s；μs, fr为固相摩擦黏度，Pa ·s。其中，固相摩擦黏度只有当固含率超过摩擦填充极限时才会考虑。
Konijn等[10]的研究表明，悬浮液体系的黏度会随着固含率的增大而显著增大。而商用Fluent软件中基于颗粒动力学理论计算得到的固相当量黏度相比于液相黏度非常小，若通过式(4)~(5)计算悬浮液黏度，在低固含率时与实际相比偏差可能并不明显，但是随着固含率的增大，造成的偏差也会随之增大。在商用Fluent软件中，搅拌功率通过扭矩与转速相乘计算得到，而扭矩由两部分组成，一部分来自于搅拌介质的压力，另一部分来自于黏性应力。因此，若将Euler模型与Mixture模型直接用于高固含率固-液两相体系的搅拌功率模拟，其准确度难以得到保证，需对其固-液两相体系黏度计算模型进行修正。
选择对Mixture模型中的固-液两相体系黏度计算模型进行修正有两点原因，一是Mixture模型将固-液两相体系作为整体进行考虑，便于实验对模拟中所需的各个参数(尤其是剪切速率)的测量；二是直接对固-液两相体系黏度计算模型进行修正较为困难，而Mixture模型中的固相当量黏度仅仅用于计算固-液两相体系的黏度(或黏性应力)，并不影响其他各项的计算，因此可以通过对固相当量黏度修正来达到修正固-液两相体系黏度计算模型的目的。通过式(5)可以得到CFD模拟中固相当量黏度与固-液两相体系黏度之间的关系如式(7)。
| $ \mu_{\mathrm{s}}=\frac{\mu_{\mathrm{m}}-\varphi_1 \mu_{\mathrm{l}}}{\varphi_{\mathrm{s}}} $ | (7) |
式中：μm为固-液两相体系的黏度，Pa ·s。需要说明的是，修正后的μs并非是实际意义上的固相当量黏度，它的修正仅仅为了得到正确的固-液两相体系黏度。
在Fluent多相流设置中固相当量黏度窗口处读取用户自定义函数(UDF)，使用DEFINE_PROPERTY宏函数作为UDF代码的主体框架。参照Fluent中的非牛顿流体黏度模型，为了防止固-液两相体系的黏度无限大或出现不合理的过小值，需要规定黏度的上下限。
2.3.3 模拟具体设置所有模拟均采用多参考系模型模拟内部动流体域(包括桨叶)的旋转，搅拌釜内的所有三维数据均采用基于压力的SIMPLE算法，湍流模型采用RNG k-ε模型，近壁面函数采用非平衡壁面函数，曳力模型选用Syamlal-obrien模型。压力的空间离散采用体积力加权方案，湍动能、动量及体积分数的空间离散均设置为QUICK，湍流耗散率的离散格式采用二阶迎风格式。时间步长初始设置为0. 001 s，等到各项残差值足够小时再改为0. 01 s用以减少计算时间。固体颗粒的最大体积分数设置为0. 61，而在设置模拟的初始条件时，在初始化后将固相以0. 55的体积分数定义Patch在流体域的最上方，这个区域对应2. 2. 2节中所描述的固-液混合层。
值得注意的是，采用Euler模型与Mixture模型计算时，固相当量黏度计算模型选用Gidaspow模型，摩擦填充极限设置为0. 6；而在采用修正后的Mixture模型计算时，除了将固相当量黏度模型替换为修正后的模型，还应将摩擦填充极限设置为大于固相最大体积分数，这是因为修正后的模型是在综合考虑了所有黏度的基础上建立的，不需要再去单独考虑摩擦黏度。
2.4 网格无关性验证为排除网格数量对计算精度的影响，分别使用145 700、277 574和428 068这3种数量的网格进行数值模拟。在釜底中心点和侧壁任一最高点之间的路径上取200个点，路径上的固相分布情况如图 2所示。
|
图 2 不同网格数下的固含率沿路径上各点分布图 Fig.2 Distribution of solid contents along the path under different grid numbers |
从图 2可以看出网格数量为277 574和428 068时，路径上的固含率分布较为相似。但是当网格数量为145 700时，固含率在路径上几个局部区域的分布与其他两种数量网格下的模拟结果差异性较大，因此，后续模拟中采用的网格数量为277 574。
3 结果与讨论 3.1 黏度测量实验结果为了更好地展现悬浮液体系黏度与固含率以及剪切速率之间的关系，并为后续的黏度模型修正提供依据，通过实验测得了不同剪切速率及不同固含率下的悬浮液黏度，如图 3所示。
|
图 3 不同固含率悬浮液黏度与剪切速率关系 Fig.3 Relationship between viscosity and shear rate of suspension under different solid contents |
从图 3可以看出，固含率为0(纯甘油)、0. 1以及0. 2时，悬浮液体系表现为牛顿流体特征，使用水平线拟合后发现其黏度分别为0. 567、0. 682、0. 999 Pa ·s。
牛顿悬浮液的黏度取决于连续相的黏度、固相的黏度、固体颗粒的形状和尺寸分布、温度以及可能存在的电荷或电场。Einstein[11]和Batchelor[12]模型表明对于硬球颗粒和中等浓度悬浮液的稀释悬浮液，悬浮液的黏度可以通过颗粒体积分数来预测，具体表达式如式(8)~(9)。
| $ \mu_{\mathrm{r}}=1+2.5 \varphi_{\mathrm{s}} $ | (8) |
| $ \mu_{\mathrm{r}}=1+2.5 \varphi_{\mathrm{s}}+C \varphi_{\mathrm{s}}^2 $ | (9) |
式中：μr为相对黏度，μr=μm/μl。
Batchelor方程在Einstein方程中增加了一个二阶项，以解释当固体颗粒体积增加时，稀释系统外颗粒之间的相互作用。根据基本假设和推导方法，常数C在4. 8~6. 2[13-14]，其中Ball等[14]发现C=5. 2对非胶体系统最准确。
因此在本研究中，对于固含率在0~0. 2内的悬浮液的黏度表达式确定为式(10)。
| $ \mu_{\mathrm{m}}=\mu_1\left(1+2.5 \varphi_{\mathrm{s}}+5.2 \varphi_{\mathrm{s}}^2\right) $ | (10) |
通过式(10)计算得到，固含率为0. 1和0. 2时，悬浮液黏度分别为0. 738 Pa ·s和0. 970 Pa ·s，与实验得到的0. 682 Pa ·s和0. 999 Pa ·s相对误差分别为8%和3%，在可接受的范围内。结果存在误差的主要原因是配置悬浮液时无法保证其完全均匀。
3.1.2 高固含率黏度关系式从图 3可以看出，高固含率的悬浮液呈现非牛顿流体的特征。影响非牛顿悬浮液黏度的因素与影响牛顿悬浮液的因素相同。此外，它还取决于剪切应力或剪切速率。Krieger和Dougherty[15]提出将高浓度悬浮液的黏度与颗粒体积分数联系起来，具体关系如式(11)。
| $ \mu_{\mathrm{r}, \gamma}=\frac{\mu_{\mathrm{m}}}{\mu_{\mathrm{l}}}=\left(1-\frac{\varphi_{\mathrm{s}}}{\varphi_{\mathrm{m}}}\right)^{-K} $ | (11) |
式中：μr, γ为恒定剪切速率的相对黏度；φm为颗粒堆积的最大体积分数；K为拟合参数。
然而，式(11)中的参数φm和K很难通过实验数据的回归得到。由于K与式(8)中φ前面的系数有关，为了便于分析，K取2. 5。当K为2. 5时，与试验数据一致，模型就失去了稳定性，因此，为了保证模型的可靠性，Liu[16]增加了一个调整参数来修改式(11)，修改后的式子如式(12)。
| $ \mu_{\mathrm{r}, \gamma}=\frac{\mu_{\mathrm{m}}}{\mu_1}=k\left(1-\frac{\varphi_{\mathrm{s}}}{\varphi_{\mathrm{m}}}\right)^{-2.5} $ | (12) |
式中：k为与粒度相关的参数。
颗粒堆积的最大体积分数φm与颗粒粒径级配、形状等因素有关，单分散体系的最大体积分数一直受到研究人员的关注，但人们对其还没有完全了解。关于这个问题已经做了大量的理论和实验工作，并且在以前的研究中提出并使用了广泛的值，大部分位于0. 55~0. 68[17]，并非严格等于随机密实堆积体积分数0. 64[18]。实验数据显示存在多组k与φm能使得拟合值与实验值取得较好的一致，一般对于微米级的球状颗粒最大堆积体积分数取0. 61[19]，因此，合理地取φm为0. 61，拟合出对应的k值为0. 86，为便于表示，定义参数ε=(1-φs/φm)-2. 5，使用剪切速率为3 s-1时的数据进行拟合，拟合具体情况如图 4所示。
|
图 4 相对黏度与ε的拟合关系图 Fig.4 Fitting relationship between relative viscosity and ε |
然而，式(12)是在某一指定剪切速率下的关系式，并不能同时反映黏度与固含率以及剪切速率这两者之间的关系。参照Liu等[20]在拟合悬浮液关系式中所作的处理，定义一个参数δ，其数值关系满足式(13)。
| $ \mu_{\mathrm{r}}=(1-\delta) \mu_{\mathrm{r}, 3} $ | (13) |
式中：μr, 3为剪切速率为3 s-1时的相对黏度，即初始相对黏度；δ为拟合参数。
为使拟合公式在固含率为0. 2~0. 3时也具有较高准确度，固含率为0. 2时的数据也参与拟合。不同固含率悬浮液的δ与剪切速率之间的拟合关系如图 5所示。
|
图 5 不同固含率悬浮液的δ与剪切速率之间的拟合关系图 Fig.5 Fitting relationship between δ and shear rate of suspensions with different solid contents |
图 5中用来拟合的关系式形式为δ=aγb+c，悬浮液拟合函数中的a、b、c都随着固含率变化，因此分别拟合它们之间的关系如图 6所示。由图 6中的拟合关系可以得到描述剪切速率与固含率之间关系的具体关系式，如式(14)。
|
图 6 参数a、b、c与固含率之间的拟合关系图 Fig.6 Fitting relationship between parameters a, b, c and volume solid content |
| $ \delta=\left(-6.27 \varphi_{\mathrm{s}}+0.92\right) \gamma^{-1.97 \varphi+0.42}-1.34 \varphi_{\mathrm{s}}{ }^{-0.33}+2.49 $ | (14) |
结合式(12)~(14)可以得到悬浮液相对黏度与固含率以及剪切速率之间的关系式如式(15)。
| $ \mu_{\mathrm{r}}=\left[\left(6.27 \varphi_{\mathrm{s}}-0.92\right) \gamma^{-1.97 \varphi+0.42}+1.34 \varphi_{\mathrm{s}}^{-0.33}-1.49\right]\left(1-\frac{\varphi_{\mathrm{s}}}{\varphi_{\mathrm{m}}}\right)^{-2.5} $ | (15) |
对比通过式(15)计算得到的相对黏度的预测值与实验值，平均绝对相对误差(MARE)为4. 17%，具体情况如图 7所示。
|
图 7 相对黏度预测值与实验值对比图 Fig.7 Comparison of predicted and experimental values of relative viscosity |
基于悬浮液黏度测量实验提出新的悬浮液体系的黏度计算模型，其目的是基于此模型修正CFD中的固-液两相体系的黏度计算模型。由于难以直接对商用Fluent软件中的固-液两相体系的黏度计算模型进行修正，因此，提出通过修正固相当量黏度进而修正固-液两相体系黏度计算模型的思路。结合式(7)、(10)、(15)可得到修正后的固相当量黏度μs的具体关系式如式(16)。
| $ \mu_{\mathrm{s}}=\left\{\begin{array}{c} \frac{1}{\mu_1}\left(3.5+5.2 \varphi_{\mathrm{s}}\right), 0<\varphi_{\mathrm{s}} \leqslant 0.2 \\ \frac{\mu_1}{\varphi_{\mathrm{s}}}\left\{\left[\left(6.27 \varphi_{\mathrm{s}}-0.92\right) \gamma^{-1.97 \varphi+0.42}+1.34 \varphi_{\mathrm{s}}^{-0.33}-1.49\right]\left(1-\frac{\varphi_{\mathrm{s}}}{\varphi_{\mathrm{m}}}\right)^{-2.5}+\varphi_{\mathrm{s}}-1\right\}, 0.2<\varphi_{\mathrm{s}} \leqslant 0.61 \end{array}\right. $ | (16) |
为探究Euler-Euler方法在模拟预测固-液搅拌功率时的具体表现，分别使用Mixture模型与Euler模型对固含率φs在0. 1~0. 4、转速N在150~300 r ·min-1的固-液搅拌功率进行模拟预测，图 8为搅拌功率的模拟结果与实验结果之间的对比图。
|
图 8 修正前搅拌功率模拟值与实验值对比图 Fig.8 Comparison of simulated and experimental values of mixing power before correction |
从图 8的对比结果可以看出，使用Euler模型和Mixture模型对固-液搅拌功率的模拟结果都存在两个问题，一是虽然固含率较低时，搅拌功率实验值与模拟值之间误差较小，但是随着固含率的增大，实验值与模拟值之间的偏差也会随之显著增大；二是搅拌功率模拟值与固含率之间呈负相关，与实验结果不相符。
产生这一情况的原因前文中已经提及过，Mixture模型与Euler模型中的固-液两相体系黏度计算模型存在不足，才使得这两种多相流模型直接应用于高固含率固-液搅拌功率的模拟预测时，模拟值与实验值的偏差较大。这一观点也可以通过功率模拟结果中的黏性应力贡献项与固含率的关系来进一步佐证，因为在CFD模拟中，功率通过搅拌的扭矩与转速来计算得到，而扭矩又是通过搅拌介质的压力以及黏性应力的作用产生的。对于压力产生的扭矩与固含率之间的关系并不好判断，但是对于黏性应力产生的扭矩与固含率之间的关系，是能够推断出的。
图 9(a)~(b)分别为使用Mixture模型与Euler模型对固-液搅拌功率的模拟值中黏性应力贡献项Pτ与固含率φs以及转速N之间的关系。从图 3的黏度测量实验结果可以看出，固-液两相体系的黏度随着固含率的增大而增大，由此可以推断出，由黏性应力引起的搅拌功率也随着固含率的增大而增大。但从图 9可以看出，黏性应力引起的搅拌功率模拟值随着固含率的增大而减小，进一步证明了Mixture模型与Euler模型中对固-液两相体系的黏度计算模型存在不足。
|
图 9 搅拌功率模拟值中的黏性应力项与固含率、转速关系图 Fig.9 Relationship between viscous stress and solid content/rotation speed in the simulated values of stirring power |
为提高对固-液搅拌功率的预测精度，以式(16)来修正Mixture模型中的固相当量黏度计算模型，进而达到修正Mixture模型中固-液两相体系黏度计算模型的目的。将进行黏度模型修正后的Mixture模型对固-液搅拌功率的模拟与实验结果进行对比，具体情况如图 10所示。
|
图 10 修正后搅拌功率模拟值与实验值对比图 Fig.10 Comparison of simulated and experimental values of mixing power after correction |
从图 10可以看出，在对固-液两相体系黏度计算模型进行修正后，搅拌功率的模拟预测精度得到很大提升，MARE仅为4. 6%。此外，搅拌功率的模拟值也与固含率之间呈正相关，与实验结果一致。
搅拌功率实验值与模拟值存在误差的主要原因是搅拌产生的部分热能在釜内没有及时散发出去，釜内介质存在一个小幅度的温度上升，而在对悬浮液黏度测量时选择的是一个固定的中间温度。
3.2.2 固相分布以及固含率测量结果进行黏度修正后的Mixture模型应用于固-液两相体系进行数值模拟的可靠性，也能进一步通过固相分布的模拟与实验结果对比来进行验证。图 11为转速为150 r ·min-1且固含率在0. 1~0. 4时，固相分布的模拟与实验结果的对比情况。
|
图 11 固相分布模拟与实验结果对比图 Fig.11 Comparison of simulation and experimental results of solid content distribution |
搅拌足够长的时间后，待釜内固相分布达到稳定，从图 11的实验图片可以较明显看出，相较于固-液混合区的浑浊，4种工况下釜底均呈现几乎无固体颗粒存在的清液区，而模拟结果也同样符合这一特征。此外，为了定量评估非清液区固相分布的实验和模拟情况，在距釜底65、100、135、170 mm的排液口处进行取样，测量出固含率具体数值后和模拟结果进行对比，对比结果如图 12所示。
|
图 12 不同高度釜壁处固含率实验与模拟对比图 Fig.12 Comparison of experiment and simulation results of solid content rate at different heights of the kettle wall |
从图 12可以看出，在非清液区模拟结果显示不同高度处的固含率数值几乎完全相同，说明固相分布几乎达到了完全均匀，这与实验的结果相近。
4 结论通过实验和模拟相结合的方法，在甘油-PMMA悬浮液体系黏度测量实验的基础上提出固-液两相体系的黏度模型，并将其与多相流Mixture模型相结合，对搅拌釜内高固含率下固-液两相体系的搅拌功率进行模拟预测，得到以下结论。
(1) 对于甘油-PMMA悬浮液体系，当固含率小于0. 2时，剪切速率对悬浮液体系黏度的作用几乎可忽略，然而当固含率大于0. 2时，甘油-PMMA悬浮液体系呈现出明显的非牛顿流体特性，剪切速率对悬浮液体系黏度的影响明显增大。
(2) 不论是Euler模型还是Mixture模型，对高固含率下搅拌功率的预测均存在较大误差，且搅拌功率模拟值与固含率之间呈现负相关。通过模拟与实验结果对比，综合分析后发现，主要原因是商用Fluent软件中固-液两相体系的黏度计算模型不适用于高固含率的情况，需要进行修正。在此基础上，提出了基于悬浮液黏度测量实验的结果进行线性回归，从而得到悬浮液体系的黏度计算关系式，再以此为依据修正商用的Fluent软件中的固-液两相体系的黏度计算模型这一思路方法。
(3) 本研究给出的固-液两相体系的黏度修正模型，应用于商用的Fluent软件中，对固-液两相体系的搅拌功率进行模拟预测，不仅在固含率低的体系(呈现为牛顿流体特征的固含率小于0. 2的甘油-PMMA悬浮液体系)中有较好的预测精度，且在高固含率(呈现为非牛顿流体特征)体系中同样有较好的预测精度。本研究中固含率为0. 1~0. 4的甘油-PMMA体系，按Mixture方法，平均绝对相对误差约为4. 6%，预测精度明显提高。
| [1] |
TAGHAVI M, ZADGHAFFARI R, MOGHADDAS J, et al. Experimental and CFD investigation of power consumption in a dual Rushton turbine stirred tank[J]. Chemical Engineering Researchand Design, 2011, 89(3): 280-290. DOI:10.1016/j.cherd.2010.07.006 |
| [2] |
董厚生, 魏化中, 舒安庆, 等. 搅拌槽内固液两相流的数值模拟及功率计算[J]. 化工装备技术, 2012, 33(1): 14-16. DONG H S, WEI H Z, SHU A Q, et al. Numerical simulation and power calculation of solid-liquid two-phase flow in stirred tank[J]. Chemical Equipment Technology, 2012, 33(1): 14-16. DOI:10.3969/j.issn.1007-7251.2012.01.004 |
| [3] |
杨红, 张远新, 王程祥, 等. 涡轮搅拌器搅拌特性数值模拟与实验研究[J]. 装备制造技术, 2012(3): 1-3. YANG H, ZHANG Y X, WANG C X, et al. Numerical simulation & experiment research on stirred characteristics of turbine agitator[J]. Equipment Manufacturing Technology, 2012(3): 1-3. DOI:10.3969/j.issn.1672-545X.2012.03.001 |
| [4] |
WU J, WANG S, GRAHAM L, et al. High solids concentration agitation for minerals process intensification[J]. AIChE Journal, 2011, 57(9): 2316-2324. DOI:10.1002/aic.12468 |
| [5] |
PAGLIANTI A, CARLETTI C, MONTANTE G. Liquid mixing time in dense solid-liquid stirred tanks[J]. Chemical Engineering & Technology, 2017, 40(5): 862-869. |
| [6] |
ZHAO H L, ZHANG T A, LIU Y, et al. Numerical simulations of solid-liquid stirred tank with an improved Intermig impeller[J]. AIPConference Proceedings, 2013, 1542(1): 1286-1289. |
| [7] |
LI X L, ZHAO H L, ZHANG Z M, et al. Numerical optimization for blades of Intermig impeller in solid-liquid stirred tank[J]. Chinese Journal of Chemical Engineering, 2021, 29: 57-66. DOI:10.1016/j.cjche.2020.08.044 |
| [8] |
GU D Y, YE M, LIU Z H. Computational fluid dynamics simulation of solid-liquid suspension characteristics in a stirred tank with punched circle package impellers[J]. International Journal of Chemical Reactor Engineering, 2020, 18(9): 20200026. DOI:10.1515/ijcre-2020-0026 |
| [9] |
李晨辰. 粉煤灰—氢氧化钠浆液槽流场的数值模拟研究[D]. 天津: 天津大学, 2014. LI C C. Numerical simulation study on the flow field in fly ash-sodium hydroxide slurry tank [D]. Tianjin: Tianjin University, 2014. |
| [10] |
KONIJN B J, SANDERINK O B J, KRUYT N P. Experimental study of the viscosity of suspensions: Effect of solid fraction, particle size and suspending liquid[J]. Powder Technology, 2014, 266: 61-69. DOI:10.1016/j.powtec.2014.05.044 |
| [11] |
EINSTEIN A. Investigations on the theory of the Brownian movement[M]. New York: Dover Publications, Inc, 1956.
|
| [12] |
BATCHELOR G K. The stress system in a suspension of force-free particles[J]. Journal of Fluid Mechanics, 1970, 41(3): 545-570. DOI:10.1017/S0022112070000745 |
| [13] |
MENDOZA C I, SANTAMARÍA-HOLEK I. The rheology of hard sphere suspensions at arbitrary volume fractions: An improved differential viscosity model[J]. TheJournal of Chemical Physics, 2009, 130(4): 044904. |
| [14] |
BALL R C, RICHMOND P. Dynamics of colloidal dispersions[J]. Physics and Chemistry of Liquids, 1980, 9(2): 99-116. DOI:10.1080/00319108008084770 |
| [15] |
KRIEGER I M, DOUGHERTY T J. A mechanism for non-Newtonian flow in suspensions of rigid spheres[J]. Transactions of the Society of Rheology, 1959, 3(1): 137-152. DOI:10.1122/1.548848 |
| [16] |
LIU D M. Particle packing and rheological property of highly-concentrated ceramic suspensions: φm determination and viscosity prediction[J]. Journal of Materials Science, 2000, 35(21): 5503-5507. DOI:10.1023/A:1004885432221 |
| [17] |
QI F Z, TANNER R I. Random close packing and relative viscosity of multimodal suspensions[J]. Rheologica Acta, 2012, 51(4): 289-302. DOI:10.1007/s00397-011-0597-3 |
| [18] |
KAMIEN R D, LIU A J. Why is random close packing reproducible?[J]. Physical Review Letters, 2007, 99(15): 155501. DOI:10.1103/PhysRevLett.99.155501 |
| [19] |
博德R B, 斯图沃特W E, 莱特富特E N. 传递现象 [M]. 戴干策, 戎顺熙, 石炎福译. 2版. 北京: 化学工业出版社, 2004. BRID R B, STEUART W E, LIGHTFOOT E N. Transport phenomena [M]. DAI G C, RONG S X, SHI T F, trans. 2nd ed. Beijing: Chemical Industry Press, 2004. |
| [20] |
LIU Y K, ZHANG Q S, LIU R T. Effect of particle size distribution and shear rate on relative viscosity of concentrated suspensions[J]. Rheologica Acta, 2021, 60(12): 763-774. DOI:10.1007/s00397-021-01301-4 |