有机废弃物通常是一种高粘度的剪切稀化非牛顿流体^[14-15]。由于高粘度的影响，其流动性与理想流体有所差别，其粘度通常可用幂律定律表示^[16]，见式(1)。

式中：η为非牛顿流体粘度，Pa·s；K为稠度系数，Pa·sⁿ；$ \dot{\gamma} $为剪切率, s⁻¹；n为流变指数。本研究的实验对象选择TS为14.3%的猪粪。根据已有研究^[17]可知，K和n分别为52.8 Pa·sⁿ和0.09。

Kompogas卧式厌氧发酵罐尺寸为36 m×9 m×12 m(工作高度9 m)。安装在搅拌式反应器侧壁面的搅拌电机驱动水平九级PBT搅拌桨。研究不同桨叶角度(30°、45°和60°)和转速对HSAD水力学条件的影响。将3种搅拌桨分别命名为叶轮 A、叶轮 B、叶轮 C ，不同工况记为Case A-30、Case B-45和Case C-60。搅拌桨具体参数有9层直径为8.1 m的叶轮、叶轮间距为3.6 m、桨叶厚度和高度分别为0.18和1.2 m。反应器结构如图1所示。

1)计算域和网格划分。本研究使用Solidworks 2016(Dassault Systemes, USA)建立几何模型，ANSYS Mesh19.0和Fluent 19.0(ANSYS Inc.，USA)划分网格和模拟计算HSAD中的水力学条件。Case A-30、Case B-45和Case C-60的网格数分别为2 554 612、2 325 531和 2 721 492。网格模型都采用非结构化网格，并在计算域中相对速度变化较大处细化网格以提高网格质量。

2)控制方程。

①连续性方程。对三维不可压流，连续性方程见式(2)。

式中：$ {\rho}_{\textit{α}} $为α相的密度，$ {{U}}_{\textit{α}} $为α相的速度，$ {\gamma }_{\alpha } $为α相的体积分数，$ {S}_{\alpha } $为质量源项。

②动量方程。α相的动量方程见式(3)。

式中：p为静压, $ \stackrel{̿}{{\mathrm{\tau }}_{\mathrm{\alpha }}} $为α相切应力张量，$ \overrightarrow{\mathrm{g}} $为重力，$ \overrightarrow{{\text{U}}_{\text{α}\text{β}}} $和$ \overrightarrow{{\text{U}}_{\text{βα}}} $为交界面速度，$ \overrightarrow{{\mathrm{R}}_{\mathrm{\alpha }\mathrm{\beta }}} $为各项之间相互作用力，定义见式(4)。

式中：$ \overrightarrow{{U}_{\beta }} $为β相速度；$ {K}_{\alpha \beta } $为交互相的动量交换系数，定义见式(5)。

式中：$ {\textit{γ}}_{\beta } $为β相的体积分数，$ {\textit{ρ}}_{\beta } $为β相的密度，f为曳力作用；$ {\tau}_{\beta} $为微粒弛豫时间，定义见式(6)。

式中：$ {d}_{\beta } $为β相气泡直径，$ {\mu }_{\beta } $为β相粘度。

③体积分数方程。α相的体积分数见式(7)。

式中：$ {\textit{γ}}_{\text{α}} $为α相的体积分数。

式中：n为相的个数。

3)数值模拟方法。由于HSAD物料的高粘度导致水力条件较差，沼气气泡的气-液相互作用可以忽略不计^[18]。假设反应器中固液两相分布均匀，将其视为“拟液相”。考虑到沼液上方的空气层，可以采用气-液两相模型预测厌氧消化内流动行为，并将液相设为主相，气相为次相。有效高度以下设为液相，其余为气相。所有壁面设为无滑移壁面。假设反应器内气液相温度恒为35 ℃，采用三维稳态模拟、欧拉-欧拉方法求解控制方程，得到速度、死区比例和搅拌功率等参数来分析HSAD的流动行为。为了合理描述不同流态下混合系统中的流动特性，在雷诺数小于10时，采用层流模型计算；当雷诺数大于10时，采用RNG k-$ {\epsilon} $模型计算。然后，指定MRF模型中的动静区域，并设置转动方向和转速大小。采用耦合式求解方案进行压力-速度耦合。采用高阶离散方案(如Second-Order Upwind、Compressive和MUSCL格式)提高准确性。

4)收敛准则。考虑到搅拌模型的复杂性，模型计算需满足2个条件视为收敛^[19]：首先，大部分求解参数残差满足$ \text{1×}{\text{10}}^{-5} $，而连续性残差和耗散率残差满足$ \text{1×}{\text{10}}^{-4} $；其次，监控截面x=0上速度平均值是否稳定，当速度均值波动小于5%时，可结束模拟并认为结果可靠。

通过系统内能量输入和输出来评价能量效率。假设厌氧发酵系统中死区内甲烷产率为0，其余区域为有效工作体积，以平均甲烷产率计。能量输入部分为电动机驱动搅拌桨耗费的能量和加热所需的能量，厌氧消化产生的甲烷视为输出能量部分，通过甲烷的热值折算为产出的能量。最终，净能量就等于输出能量减去输入能量。计算方法见式(9)~式(13)^[20]。

式中：$ E_{\text{mixing}} $是混合能量，GJ·d⁻¹；N是搅拌桨转速，r·min⁻¹；P是总力矩，N·m；$ E_{\text{heating}} $是加热能量，GJ·d⁻¹；ρ 是水的密度，kg·m⁻³；Q是流量，t·d⁻¹；C是水的比热，J·kg⁻¹·℃⁻¹；T_d是消化温度，℃；T_a是环境温度，℃；k是传热系数；V_工作为工作体积(有效高度区域体积)，m³；DZR为死区体积比例；$ E_{\text{output}} $是从甲烷回收的能量，GJ·d⁻¹；$ P_{\text{methane}} $是甲烷产率，m³ CH₄·kg VS⁻¹；q是甲烷的热值，J·m⁻³；OLR是有机负载率，kg VS·m⁻³·d⁻¹；V是有效工作体积，m³；$\eta $是能量转化率，%；$ E_{\text{net}} $是净能量，GJ·d⁻¹.

有研究将小试反应器的能量评估方法推广到大试规模反应器下^[3]。在前期研究中，已经使用工作容积50 L的反应器开展实验，故本研究也使用该方法进行能量评估。有关运行参数(如OLR、甲烷产率、消化温度等)的详细描述参考^[20]。

前期研究^[7]发现，CFD结果可以与PIV测量结果进行比较以验证其是否可信。3.00%的锂藻土悬浮液($ \text{η}\text{=}\text{K}{\dot{\text{γ}}}^{\text{n-1}} $，其中k=36.4 Pa·sⁿ、n=0.05)和猪粪的流变特性和结构特性相类似，被用来作为PIV的模拟工作流体^[17]。PIV测试系统与前期研究^[7]一致：在透明圆柱反应器中装置3.00%的锂藻土悬浮液，并采用双折叶搅拌桨进行连续搅拌，圆柱反应器直径为0.1 m、高度为0.2 m、有效高度为0.12 m。将上述装置置于装有水的透明方形容器内，由于3.00%的锂藻土悬浮液与水的折射率相同(均为1.33)，可以抵消光学折射产生的实验误差。CFD模拟结果可以和PIV测量结果比较来确定模拟结果是否可信。如果模拟结果被验证合理，那么其他的物理参数(如流场、死区等)也认为是可信的。

非牛顿流体机械搅拌的雷诺数(Reynolds number, Re)定义见式(14)^[21]。

式中：N为搅拌桨转速，s⁻¹；D为桨叶直径，m；$ \textit{ρ} $为厌氧消化物料的密度，kg∙m⁻³；K为稠度系数，Pa·sⁿ；n为流变指数。

搅拌功率可以由力矩计算得到，定义见式(15)^[22-23]。

式中：P为搅拌功率，W；N为搅拌桨转速，s⁻¹；T为总力矩，N·m。

转速为150 r·min⁻¹下的CFD模拟和PIV测量结果对比如图2所示，2者的流型和速度涡分布具有一致性。在CFD模拟结果中，150 r·min⁻¹下桨叶边缘处的最大速度为0.10 m·s⁻¹，而PIV测量结果显示在该区域的流速为0.09 m·s⁻¹。CFD计算速度与PIV测量值之间的差异可能是由于高速摄像存在一些光学像元和噪声^[7]，而CFD结果是一种理想化的数值模拟方法。CFD和PIV结果之间差别在10%以下可知，本研究建立的CFD模型是可靠的，能够合理反映整体的水力学条件。

1)不同桨叶角度对径向混合的影响。Case A-30、Case B-45、Case C-60在不同转速下的径向流场(z = −3.6 m)如图3所示。与桨叶尖端区域相比，其他区域在0.5 r·min⁻¹转速下的流场流速较低。由此可以推断，消化物料仅在桨叶边缘处有效混合。

当转速达到1 r·min⁻¹时，有效混合区呈等边三角形。随着桨叶角度的增大，桨叶间的混合沿径向逐渐向外扩散，进一步诱导了更多的液体流动。这说明，桨叶角的增大促进了径向混合。也就是说，在相同转速下，桨叶角的增大导致径向剪应力增大。根据猪粪的流变特性，伴随剪切应力的增大，猪粪粘度降低，进而可以提高猪粪的可流动性。其次，值得注意的是，两级叶轮之间的区域在小桨叶角时(叶轮A、叶轮B)是接近静止的，而在大桨叶角时(叶轮C)则是充分混合的。

当转速大于3 r·min⁻¹时，有效混合区域由近似三角形的形状向近似圆形的旋转形状演变。这表明，转速的增加也可以增强混合。随着转速从3 r·min⁻¹增加到7 r·min⁻¹，近似圆的直径从8 m增加到9 m。即，混合首先在搅拌桨周围产生，然后随着转速的增加，混合区域逐渐扩大到壁面，从而使反应器内的物料能够有效混合。例如，当转速为3和5 r·min⁻¹时，在x = 0~3 m、y = 0~3 m和x = 6~9 m、y = 0~3 m(本节使用的x和y坐标是基于图3)的范围内混合强度明显不足，而转速提高到7 r·min⁻¹，混合不足的区域缩小到x = 0~1 m、y = 0~2 m和x = 8~9m、y = 0~1m。

但桨叶角度对径向混合的影响与转速相比并不明显，特别是当叶轮转速大于3 r·min⁻¹时，叶片角对径向混合的改善效果非常有限。另外，在不同的工况下，这些反应器的底部角落均不能有效混合。其可能的原因是：桨叶尖端周围流体受到最大的剪切力，但远离搅拌桨的反应器底部受到的剪切力最小，故导致该区域的混合液粘度较高，不能有效混合。

2)不同桨叶角度对轴向混合的影响。Case A-30、Case B-45、Case C-60在不同转速下的轴向流场(x=0 m)如图4所示。可以看出，随着转速的增加，反应器内的流速逐渐增大，搅拌桨的扰动范围逐渐扩大，这与已有研究的观察结果一致^[7]。在Case A-30, Case B-45和Case C-60中，消化池的流场模拟显示低转速下(0.5 r·min⁻¹)存在一个很大的未混合区。这部分区域内速度接近零且粘度接近零剪切粘度，即形成死区^[24]。综上所述，转速为0.5 r·min⁻¹工况下反应器轴向混合不均匀，故物料无法有效混合。

而当转速增加到1 r·min⁻¹时，配置叶轮A、叶轮B、叶轮C的反应器流场有明显差异。不同于0.5 r·min⁻¹下的流场结果，叶轮两级间的停滞区开始发生混合，即轴向混合得到改善。另外，Case B-45和Case C-60的混合区域沿轴向扩大，并且在Case A-30的1 r·min⁻¹工况下，混合区域扩大到叶轮各级之间的区域。这一现象不同于已有研究，即叶片角的增大会促进径向混合(如3.2.1所述)，而对轴向混合的影响很小。这可能是由于桨叶角度的扩大增大了物料受到的剪切力，提高了消化液的可流动性，即提高了流动速度，从而增强了轴向混合。综上所述，叶片角对轴向和径向混合都有较大的增强作用。

同样，当转速达到3、5、7 r·min⁻¹时，叶轮C的轴向混合表现出更充分的流动性。在转速从3 r·min⁻¹增加到5 r·min⁻¹，充分混合区域为z = 3~33 m、y = 0~9 m(本节使用的x和y坐标是基于图4)，扩大到x = 2~35 m、y = 0~9 m。其次，随着转速增加到7 r·min⁻¹，叶片角的增加也会增强轴向混合，但作用范围很有限。此外，从图4中可以看出，相邻叶轮接近底部区域总是难以有效混合，如图3(c)所示呈现“三角区”。这可以解释为：左边叶轮驱动液体向右上方运动，而右叶轮驱动液体向左上方向运动，使得“三角区”的流体受到的作用力抵消，产生停滞现象。

死区可以用来评价流动模式和混合效率。定义反应器流场速度低于0.001 m·s⁻¹的部分标记为死体积^[25]。假设死区为不发生流动的区域^[25]，该区域的对流传质效率接近0。除死区外，其余区域划分为3个速度区，分别为低速区(Low velocity zone, LVZ)(0.001 ~ 0.01 m·s⁻¹)、中速区(Medium velocity zone, MVZ)(0.01 ~ 0.1 m·s⁻¹)和高速区(High velocity zone, HVZ)(>0.1 m·s⁻¹)^[26]。

转速和桨叶角度对流速分布的影响如图5所示。可以看出，在3种构型的低转速(0.5和1 r·min⁻¹)下，MVZ所占体积分数最高。但随着转速的不断增加，该比例逐渐降低。然后，随着转速在3 r·min⁻¹以上继续增加，流速在0.1~ 1 m·s⁻¹区域占主导地位。此外，当转速高于3 r·min⁻¹时，LVZ的体积分数很低(小于14%)。另一方面，在0.5和1 r·min⁻¹的低转速下，HVZ所占体积分数相当低，但随着转速增加到3 r·min⁻¹，HVZ逐渐占主导地位。另外，在Case B-45和Case C-60中，转速到达7 r·min⁻¹条件下，HVZ中流速低于1 m·s⁻¹的区域开始减少，但流速高于1 m·s⁻¹的区域迅速增加。

从图5可以看出，桨叶角度对流速分布也有很大影响。特别是在低转速下，桨叶角度的影响比高转速下更为明显。例如，死区和LVZ比例随着叶片角度的扩大显著减少，而MVZ和HVZ在0.5 r·min⁻¹转速下增加。另外，MVZ的体积分数顺序为：45° > 60° > 30°。HVZ的体积分数顺序为：60° > 45° > 30°。随着转速增加到1 r·min⁻¹以上，死区和LVZ呈现出相同的趋势。而MVZ则随着角度的增加而减小，随后出现速度大于1 m·s⁻¹的HVZ。在3种情况下，可以明显看出，死区和LVZ随着转速的增加而逐渐减小，HVZ中速度高于1 m·s⁻¹的区域随着转速的增加而逐渐增大。

Case A-30、Case B-45和Case C-60的流速空间分布如图6所示。可以看出，死区只存在于靠近反应器壁的区域，并且死区保持在0.28%~5.65%的比率，这与其他研究^[4]一致。此外，可以清楚地看出，在0.5和1 r·min⁻¹转速下，角度的增加增强了反应器上部区域的气体流动，但对流场的改善作用不大。当转速超过1 r·min⁻¹时，桨叶的扩大只增大了叶轮周围区域的流速，对死区体积减小的影响可以忽略。因此，在评价大试HSAD反应器的流态时，不仅需要考虑死区体积，还需要考虑LVZ、MVZ和HVZ的比例。

如图7(a)所示，搅拌桨消耗的能量随着转速的增加而增加，当转速从0.5 r·min⁻¹增加到3 r·min⁻¹时，不同桨叶角度消耗的能量差距相对很小。然而，随着转速继续增加到7 r·min⁻¹时，混合消耗的能量差异愈发显著。不同角度叶轮所消耗的混合能量高低顺序为：30° < 45° < 60°。

图7(b)为叶片角度和转速对死区的影响。由图可知，死区比例随转速的增加而减小，但减小的幅度很小，转速从0.25 r·min⁻¹增加到7 r·min⁻¹时，死区比例仅减小不到5.5%。相较于30°和45°，60°可以在低转速下达到相同的死区比例。在0.25和1 r·min⁻¹的低转速下，桨叶角度的扩大与死区的减小呈现一致性；而在高转速下，这种影响可以忽略。这说明，从混合能量的经济性和低死区考虑，均不需要高转速。

输出能量指从沼气中回收的能量，与反应器内的有效混合区域密切相关。图7(c)为不同工况下的输出能量情况。输出能量随着转速的增加而增加。其中，在0.25和1 r·min⁻¹的低转速下，桨叶角度对输出能量的影响相当明显，而在高转速下可以忽略。这种趋势与混合能量的趋势差异较大。

不同工况下产生的净能量如图7(d)所示。在不同桨叶角度下净能量变化趋势相似。由图中可以观察到，在0.25 r·min⁻¹转速下配备60°叶轮的混合反应器在3个反应器中表现出最高的能量性能，其净能量为66.823 GJ·d⁻¹。相比之下，装有30°和45°叶轮的搅拌反应器产生净能量分别为64.727和65.349 GJ·d⁻¹。在0.25和1 r·min⁻¹的低转速下，净能量处于稳定水平。然后随着转速继续增大，净能量急剧下降。这很大可能是由于提高转速增加了混合能量，但其增强流动性的作用非常有限。

综上可知，从净输出能量的角度来看，在转速为0.5和1 r·min⁻¹时，配备60°叶轮的混合系统获得了最高的净能量。从工程应用方面来看，在转速为0.5 r·min⁻¹时，电机效率没有得到充分利用。因此，推荐采用60°叶轮，转速为1 r·min⁻¹的HSAD反应器。

对最优工况(叶片角为60°，1 r·min⁻¹)下的流速进行定量评价。叶轮叶尖、中部、搅拌轴轴向混合流场如图8(a)~ 8(c)所示，相对位置如图9所示。可以看出，沿搅拌轴对称位置的流速不是对称的，这可能是由于雷诺平均模型的影响。可以看出，不同叶片的轴向速度均表现出从叶轮周围到两级叶轮之间区域逐渐减小的趋势。这主要是由于叶片周围的流体受到较大的离心力，推动流体运动。这种离心力随着轴向距离的增加而逐渐减弱。在轴向中间位置，污泥仅受到轻微扰动，流速维持在0.01 m·s⁻¹左右。靠近搅拌轴区域的相对流速小于0.004 m·s⁻¹。由此可以推断，搅拌轴只对物料的流动提供摩擦力，不能有效地诱导物料的流动。

不同径向位置流场分布如图8(d)所示。α-30°、α-60°、α-90°径向线流速呈现先增大后减小的趋势。但是，速度峰值出现的位置不同，α-30°径向线在x = 1 m处出现峰值，α-60°径向线在x = 1.7 m处出现峰值，α-90°径向线在x = 2.5 m处出现峰值。因此，可以得出结论，流场在轴向不同位置的分布是不均匀的。

1)基于CFD实现HSAD流场可视化发现，桨叶角度从30°扩大到60°可以明显改善径向混合和轴向混合。

2)在较低转速下(0.5 r·min⁻¹或1 r·min⁻¹)扩大桨叶角度可以显著减小死区。建议采用低转速和大桨叶角度。

3)针对Kompogas卧式厌氧发酵罐，推荐使用桨叶角度60°的斜叶式搅拌桨在1 r·min⁻¹转速下以获得最优的净能量输出0.023 GJ·d⁻¹·m⁻³。