大型发酵机多层搅拌流固耦合仿真分析

徐帅龙，叶建华，徐欢，唐辉，曾寿金，韦铁平

(1.福建省智能加工技术及装备重点实验室，福建 福州 350118；2.福建工程学院 机械与汽车工程学院，福建 福州 350118)

大型封闭式高温好氧发酵处理机正成为畜禽废弃物无害处理、资源再利用[1-3]的主流设备。当前，由于大规模集约化养殖的要求和土地资源的日益紧张，发酵设备需求正向100 m3以上的大型容积方向发展。搅拌刀具是大型封闭式发酵机的核心部件。在大型塔式密闭环境下，传统的单层刀具难以满足半径大、堆深高的高温发酵搅拌要求。多层搅拌刀具具备分层搅拌的能力，能满足大型密闭环境下的搅拌需求，但是容易出现负载过大导致失效或者搅拌混合效果不佳等问题。因此有必要对搅拌刀具进行流场和流固耦合仿真分析，研究搅拌刀具的搅拌效果和承载能力，进而保障整机性能。

在刀具的搅拌混合效果和受力变形方面，学者们从不同的角度进行了相关研究。钱皎龙等[4]通过改变刀具离底高度开展研究，认为提高搅拌刀具离底高度有利于减小流场内搅拌刀具的死区、提高其混合效果。罗松等[5]研究发现把四直叶圆盘涡轮桨固定在任意位置时，将六折叶开启涡轮桨片安装在二折叶桨的上方时整体结构的混合时间更短。刘敏珊等[6]通过对斜叶涡轮桨的3个不同角度对比搅拌刀具不同倾斜角度对搅拌罐内的混合响应情况。唐睿等[7]对卧式双轴发酵机更换不同的材料属性对比观察相同工况下的搅拌装置转矩、功率、以及所受应力和应变情况。巩亚东等[8]从搅拌物料的黏度属性角度出发，对比模拟其对搅拌刀具受力、变形的影响情况。在流固耦合作用下，刀具转速也是影响混合效果以及刀具受力、变形不可忽视的因素。针对非牛顿流体，方玉建等[9]对机械搅拌罐内非牛顿流体内流特性进行研究，发现提高转速，搅拌槽内涡的分布和湍流动能分布变大，轴向速度和混合时间也相应缩小。

综上，本研究基于畜禽养殖废弃物的牛顿流体特性，对流固耦合作用下的刀具转速与混合效果及刀具受力、变形的影响关系开展研究，为搅拌发酵设备的标准制定以及工艺参数选取提供参考和借鉴。

1 物理模型

1.1 搅拌槽及刀具结构

根据规模化养殖场禽畜污物的处理能力需求，设计发酵设备的容积为116 m3，采用相似理论对原尺寸进行三维软件建模，发酵罐的结构示意图如图1所示。

图1 发酵罐结构示意图(单位：mm)

搅拌轴由短轴和长轴通过连接法兰相连而成，由于轴长超过4 m，故用空心轴[10]代替实心轴以提高搅拌系统的临界转速，节约了大量钢材。为获得更好的搅拌效果并降低功耗，将搅拌刀具设计成多层结构。其中，最上面一层采用两斜叶开启涡轮式结构，下面两层为四斜叶开启涡轮式结构，叶片尺寸为280 mm×35 mm×5 mm，安装倾斜角度为45°。

1.2 运行工况及边界条件

搅拌过程中，搅拌轴的顶部和底部分别采用调心和推力轴承固定，对搅拌轴底端设置固定支撑和圆柱支撑约束。由于搅拌罐壁及罐顶在投入搅拌物料后采用封闭处理且罐顶受搅拌扰动较小，将罐壁边界条件设置成Wall，罐顶边界条件设置为Symmetry[11]，载荷示意图如图2。发酵机一次发酵处理的物料含水率较高(大致为60%)，物料黏度较低，属于低粘度流体，故其黏度值取5 Pa·s。发酵设备任务指标为日处理量大于14 m3且处理周期小于等于6 d。由于搅拌黏度较低且期望任务量较大，为了取得较好的搅拌效果，取不同的较大转速对搅拌罐进行仿真分析。

图2 载荷/约束示意图

2 数值模拟

2.1 流场仿真

2.1.1 搅拌槽内介质及流动形式

用于数值模拟的搅拌罐内的搅拌物料为粪便畜污与水的混合物，为了方便罐内物料与好氧菌的混合发酵，物料按一定的比例混合，其混合物属性接近牛顿流体的流变特性，流动雷诺数Re通过式(1)计算。

(1)

式中，ρ为介质密度，经调定为800 kg/m3；μ为动力黏度，取5 Pa·s；D为搅拌叶片的直径，取610 mm；N为搅拌叶片的转速，对N选取100、200、300、350 r/min 4个较大的转速进行数值模拟。4种转速下所得到的雷诺数值Re分别为99、198、297、347，均偏小，因此认为搅拌罐内的流体流动皆属于层流状态。

2.1.2 网格划分及材料定义

利用高级流体仿真前处理工具ANSYS-Fluent meshing对流体区域划分网格。由于桨叶附近区域几何结构复杂、介质对流强烈，因此转子区需要对网格进行局部加密，搅拌桨叶动璧面进行二次加密。网格节点共计1 895 847个、网格单元有701 511个。在ANSYS的Workbench模块中对静力学和动力学进行连接搭建，将搅拌桨模型导入到Workbench的Static Structural模块中，在Static Structural中对各部件进行材料属性的赋予。其中搅拌轴采用Q235材料，密度为7 830 kg/m3、泊松比为0.274、弹性模量为210 GPa；搅拌桨采用Q345材料，密度为7 850 kg/m3、泊松比为0.280、弹性模量为206 GPa。然后对搅拌轴和搅拌叶片采用与流场网格划分相同的尺寸进行四面体非结构化网格划分，最终确定的网格单元总数为39 580个。

2.1.3 流场求解设定

数值模拟采用多重参考系法(MRF)，首先将搅拌罐内的流体区域分为两个部分，分别为旋转的转子区和静止的定子区，再对转子区和定子区使用不同的坐标系求解。转子区和定子区之间通过内部交界面进行数据的交换；搅拌罐内流场的数值计算选用分离隐式求解器；粘性模型选择层流状态；压力与速度的耦合采用半隐算法；压力离散项采用标准格式；动量离散选择一阶；将流场进行标准初始化；设置迭代步进行迭代计算，直至收敛。

在搅拌罐的一处加入示踪剂，当采用示踪剂浓度法跟踪模拟搅拌物料在不同转速下的混合效果计算其混合时间时，需要用到组分质量守恒定律，如式(2)。

(2)

式中，ρcs即该组分的质量浓度；cs为组分s的体积浓度；Ss为系统内部单位时间内单位体积通过化学反应产生的组分s的质量；Ds为组分s的扩散系数。

要进行非稳态的浓度场计算，首先需要得到稳态下的流场数据结果，再将其作为基础条件，然后打开组分模型、激活组分传输项。在结果控制中关闭能量和流动方程，保证搅拌罐内只进行单纯的示踪剂浓度求解过程，设定浓度残差为10-4时达到收敛，时间步尺寸设置为0.02 s，最大迭代步设置为20。

3 仿真结果分析

3.1 流场分析

混合搅拌罐内物料所需的能量主要由叶片提供。通过搅拌叶片的剪切作用使得周围的搅拌物料获得能量，并以一定的规律搅动物料向罐内其他区域做扩散运动，从而达到混合的效果。取中间转速N=200 r/min为例，对搅拌叶片作用下的流场进行分析。图3为N=200 r/min时的流场速度云图和矢量图。

图3 N=200 r/min速度云图和矢量图

由图3(a)可见，最大速度出现在叶片端部，桨叶附近及双层桨间速度较大，转轴下方、罐壁和四周罐底处易形成死区。由图3(b)可见，上、中、下3层桨叶在叶片末端产生水平射流射向搅拌罐壁，从而带动附近的流体向四周运动。流体与搅拌罐的四周发生碰撞后分别向上和向下分流，然后转变为流向转轴的径向流动，最后形成回流返回搅拌叶片，即形成涡流。其中上桨由于距离中桨较近主要形成轴向流；中桨和下桨各形成4个对称的大涡流，可以起到一个将罐底的物料充分搅拌的效果，主要表现形式为径向流。图3(c)、3(d)是下层桨水平面即Z=310 mm的速度云图和矢量图，通过观察发现水平面的速度呈一定规律分布，在半径为300 mm处出现最大值。速度集中分布在搅拌桨的尾翼附近，沿着其转动方向在一定范围内逐渐减小，在趋势上由桨叶末端向转轴中心和罐壁方向逐渐减小。

因此，搅拌罐内最大流速位于底层桨端部附近区域，最大速度值为6.41 m/s；最小速度位于搅拌罐内壁面和搅拌桨底部区域以及搅拌轴附近位置；N=200 r/min时，能实现全罐搅拌。

3.2 不同转速对混合效果的影响

为了研究转速对混合效果的影响，在搅拌罐内不同位置处设置监测点对示踪剂的浓度进行实时监控来反应搅拌桨的混合时间。混合时间[12]T0被定义为自搅拌过程开始到搅拌罐内物料的物理或化学特性参数不存在明显差异的时长。国际上一般采用95%的规则来定义混合时间，即搅拌罐内一个或多个监测点达到最终稳定浓度值的±5%时所使用的时间可认为是搅拌混合时间。由于本研究分析各搅拌桨相对混合时间，对全罐混合要求不高，仅选取介于搅拌桨附近位置的一个投料点和两个监测点对搅拌罐进行监控，具体位置如图4所示。

图4 投料点和监测点位置分布

在投料点加入示踪剂，示踪剂随着搅拌物料的流动进行罐内扩散，当搅拌罐内所有位置的示踪剂浓度都几乎相等时，即可认为混合完成。由于桨叶末端处速度最大，扩散速度最快，故选择在最常用的加料点T(桨叶末端处)加料。在加料点T加料，观察不同转速下，搅拌罐内X=0铅垂面示踪剂浓度随时间变化的过程如图5所示。

由图5可见搅拌罐内以任何转速运转都会出现一个“隔断”现象，示踪剂会集中在二层桨叶附近并出现聚集现象。随着示踪剂在上半层逐渐稳定和上下层浓度差的形成，示踪剂被二层桨径向泵送到罐壁，发生撞击之后逐渐往下扩散，这与图3(b)中物料的流动特性保持高度一致。示踪剂的扩散速率会受到介质运动影响。随着混合过程的进行，上下两层的浓度差越来越小，扩散动力逐渐减小，直至两层浓度近乎相等，混合过程即可认为完成。

由于搅拌叶片的不同功用特征，示踪剂扩散特性也受到影响，每个监测点所测浓度响应曲线也是大不相同。如图6所示监测点1位于三层桨末端，相较于监测点2，点1距离投料点T更近，并且位于“隔断”上方，因此监测点1的混合时间整体上要快于监测点2。观察图6(a)可以发现转速为100 r/min时的混合曲线波动幅度最大，这是因为示踪剂的浓度在初始状态是单位1且是一个不断降低的趋势，其在相同时间内扩散速率最慢；转速300 r/min和350 r/min时的扩散速率明显快于200 r/min以及100 r/min,说明在一定范围内转速是影响混合效果的重要因素。对比图5可以发现高转速下示踪剂扩散的范围更大，示踪剂在数值上降低的更快。由于图5、图6在结果上的相互印证，可以认为相同时间内高转速情况下叶片能够在发酵罐掠过更大的范围，带动更多的物料进行混合。但在此种结构尺寸搅拌系统设计下，在转速达到300 r/min时，继续增大转速对混合效果的增益不大；转速在100 r/min以下时，搅拌混合效率不高。因此在进行发酵工作时，可以根据发酵进度的需求在100～300 r/min范围内进行调整。

图5 不同转速下示踪剂浓度扩散分布图

图6 不同转速下示踪剂浓度曲线

3.3 不同转速下搅拌装置变形情况

发酵罐内搅拌物料在不同转速下，搅拌桨每层叶片最大变形量如图7所示。分析图7可知在搅拌装置总体变形中，底层桨叶末端变形位移量最大；中层桨叶末端位移变形量也比较大；上层搅拌桨叶末端变形量最小。每层搅拌桨的最大变形量均出现在桨叶末端，随着靠近搅拌轴变形量均匀减小。由图7可见，随着转速的增加各层桨的变形量增速趋势大致相同，其中转速为200 r/min时，三层桨的最大变形量都在5 mm左右；当转速达到300 r/min时，三层桨的最大变形量达到了12 mm左右，较200 r/min变形量变化较大，说明转速200 r/min是搅拌装置的一个关键转速点。在考虑功耗不追求搅拌效率的情况下，可以将整机转速调低至200 r/min左右，其变形量相对较小，不会影响发酵机的正常工作运转。由于每层搅拌桨叶的扇数都是偶数设计，因此各层桨叶的不同搅拌桨变形趋势相同，没有出现变形不均或者某处突变的现象。

图7 不同转速下搅拌桨变形折线图

3.4 不同转速下搅拌装置受力情况

忽略搅拌装置在启动时搅拌叶片受到的瞬间冲击作用，不同转速下，搅拌装置各层桨叶所受到的等效应力分布如图8所示。搅拌装置应力分布主要集中在桨叶上，随着转速的增大，每一层桨的最大应力值都随之增大且上层桨根部所受应力值始终最大，这与上层桨所受单位力矩最大相一致。其他层桨叶所受应力最大处同样出现在桨叶根部位置，应力值沿着桨叶向外逐渐减小。

从图8发现在速度小于300 r/min时，二层桨所受应力值稍大于一层桨，随着速度的逐渐增大，中层桨的最大应力值小于底层桨。但是从总体趋势上来看，每一层桨的最大应力值都相差不大，且都小于搅拌桨材料Q345的屈服强度345 MPa，很好地避免了由于转速增大而造成某一层搅拌桨折断的问题。当速度达到350 r/min时，搅拌桨与焊接处的最大应力已经达到了261 MPa,取安全系数为1.5，则应力值为391.5 MPa，超过了搅拌桨使用的材料Q345的屈服强度。故以该速度长时间运转会存在搅拌桨折断的风险。研究结果还验证了搅拌装置转速为320 r/min时的最大应力值，发现其最大值为235 MPa，取安全系数为1.5，则应力值为352 MPa，略大于搅拌桨材料Q345的屈服强度。

图8 不同转速下搅拌桨等效应力柱状图

总之，在试验工况下，搅拌装置的极限转速为320 r/min，长期超过该转速会造成搅拌桨材料的不可逆变形甚至折断现象。因此基于搅拌稳定性上，该系统转速应调至200 r/min以下较为适宜。

4 结论

本研究使用ANSYS软件中的Fluent和Workbench模块对多层桨搅拌器进行了单向流固耦合数值模拟分析，分别得到了流场结构下结构场的数值结果以及转速对混合效果的影响。具体结论为：

1)通过对搅拌罐内的流场情况分析得知，搅拌罐内最大流速位于底层搅拌桨端部左右附近区域；最小速度位于搅拌罐内壁面和搅拌桨底部区域以及搅拌轴附近；当N=200 r/min时实现全罐搅拌。

2)在搅拌效果上，发酵罐内物料混合效果随着搅拌转速的增加提升比较明显。当转速达到300 r/min时，继续增大转速对混合效果的增益不大；转速在100 r/min以下时，搅拌混合效率不高。所以综合考虑，搅拌转速在100～300 r/min时，混合效果最佳。

3)在正常工况下，搅拌装置底层桨变形最大；各层桨变形最大处均在桨叶末端；每一层桨应力最大处均发生在彼此相互连接的轴段处，且随着搅拌转速的增加，应力应变都相应增大。当转速超过200 r/min时，搅拌叶片的最大变形量突然增大，最大应力增长也较快，但最大受力点和最大变形处不随转速的变化出现改变。因此可在应力集中的位置进行强化处理，防止搅拌叶片过载情况下失效。

综上可知，在本研究的搅拌装置作用工况下，搅拌转速设置为200 r/min左右时各项指标均能满足要求，系统比较稳定、能效最高。本研究结果能为进一步优化搅拌结构以及转速的设定提供一定的参考。