利用CFD多相流模型优化设计沼气料液搅拌流场的方法研究

黄如一， 黄正昕， 冉 毅， 熊 霞， 梅自力， 王 军， 孔垂雪

(1. 农业部沼气科学研究所 农业部农村可再生能源开发利用重点实验室， 成都 610041； 2. 四川大学 建筑与环境学院， 成都 610065； 3.上海交通大学 国际与公共事务学院， 上海 200240； 4. 四川省农村能源办公室， 成都 610041)

沼气工程通过消耗秸秆、畜禽粪便、生活污水等“三废”污染物，生产高品质清洁能源，兼具强大的环保效益、能源效益和减排效益，是解决农村能源短缺和环境污染的实用节能减排技术，无论在发达国家还是发展中国家，普遍被视为最具希望的可再生能源方式[1-3]。但目前沼气工程大多采用自然进出料和自然布水发酵，极少人工干预，亦缺乏有效的干预方法，所以普遍存在流动性差、产气率低的缺陷[4]。解决这一问题的路径是采取人工搅拌，打破静态发酵，优化流态[5]，现代沼气工程普遍都要添加搅拌工艺[6]。搅拌的具体形式主要分叶轮机械搅拌、水力搅拌、气动搅拌3种，Karim[7]通过实验证明水力搅拌是其中最优的一种。然而，由于沼气发酵装置的密闭性，料液搅拌中的流动状态难以监测，其流场形态亦是不可见的，从业人员往往只能在没有掌握流场形态的情况下，随意添加一些搅拌措施，但实际上并不适用于其罐体形状或原料特性[8]，不但不能有效提升发酵效率，甚至反而消耗更多能源[9]。

进入21世纪以来，沼气行业开始利用数值模拟方法，定量描述沼气料液搅拌时的流场形态，作为优化设计搅拌措施的依据，取得了长足的进步[10]。计算流体动力学(Computational Fluid Dynamics，CFD)软件是当前最常用的数值模拟工具，业界利用CFD方法极大地提升了沼气工程设计水平，并从实验和实际工程两个途径切实证明了CFD方法对沼气料液搅拌流场模拟计算的准确性[11]。冯琳[12]通过优化设计一种底部进水、高位分散式压力出口罐型的过程，阐述了利用CFD方法排除流场缺陷，优化设计罐型的方法，但仅限于较简单的单相流方法。事实上，沼气料液大多是复杂的多种发酵原料联合发酵[13]，并至少包含固、液两种相态，所以利用CFD的多相流模型模拟计算多种相态沼气发酵原料的搅拌过程至关重要。本文以一个底面直径0.8 m，高1.5 m的圆筒形沼气发酵罐为例，定量描述水力搅拌不同进水流速的搅拌效果，尤其是在料液呈明显多相流特征的发酵后期，得到最佳的水力搅拌速度的过程，阐述CFD多相流模型在沼气料液搅拌设计中的应用方法。

1 模拟对象和数学模型

1.1 模拟分析对象与条件

为分析水力搅拌对沼气发酵装置内部流场的具体影响，构建一个循环系统，通过泵提供动力，从底部进水，通过4个分散式出口出水，形成循环搅拌。罐体底面直径0.8 m，高度1.5 m，发酵料液有效容积为0.6 m3，储气间容积为0.15 m3。出水管高度距罐底面1.1 m，共4个。如图1所示。

1.循环管； 2.液柱； 3.储气间； 4.导气管； 5.出水口； 6.进水口； 7.水泵图1 分散式压力出口设计图

1.2 数学模型

1.2.1 以水为介质的基本流场数学模型

沼气发酵料液前期为低浓度混合溶液，流动性近似于水，所以用液态水作为介质，模拟其流场，可以作为该沼气发酵装置中流体的基本流场。

不考虑时间变化的影响，其流体流动连续性方程如下：

(1)

其动量方程如下：

(2)

(3)

(4)

式中：▽为哈密顿微分算子：

(5)

1.2.2 多相流形态采用欧拉模型

发酵后期，秸秆逐渐转变为高浓度粘性流体，呈假塑性流体相态[14]，其搅动相当于秸秆溶液与水的液-液两相流，且浓度较高。根据王军[15]给出的沼气料液搅拌工况适用模拟算法模型表，该阶段工况第二相体积分数远远超出1%，所以不宜采用拉格朗日算法下的离散相模型，而应该在欧拉算法下的3种模型中选择。考虑到该阶段秸秆原料已经成为假塑性流体而不再是颗粒形态，所以排除体积分数模型。考虑到假塑性流体(第二相)与水(第一相)的相间拽力较大，且两相之间分离不明显，所以排除混合模型，最终决定采用欧拉模型。假设两种流体相分别为p和q，其q相的流动连续性方程为如下：

(6)

其中q相的动量守恒方程：

(7)

(8)

q相的体积分数连续性方程如下：

(9)

1.3 数值方法

沼气发酵罐主要由发酵料液区、储气区两个主要部分组成，水力搅拌所影响的主要是料液部分，所以将液柱作为模拟计算对象，模型简化为1.2 m高的液柱，底部有一个直径0.05 m的圆形孔洞，1.1 m高处均匀分布4个直径0.1 m的圆形孔洞。采用控制容积法划分网格，网格类型选用非结构化四面体网格，共生成3334851个网格，如图2所示。

图2 液柱模型的网格示意图

2 不同工况的模拟结果对比

2.1 以水为介质的模拟

2.1.1 进水速度为0.3 m·s-1时的模拟结果

图3是进水流速为0.3 m·s-1时的模拟结果，从该图可见，流速分布相对均匀，且较多动能分布在高位。发酵初期较多物料浮于上方，所以动能集中在上方有利于克服重力，使其混合均匀，提升发酵效率[7]。但该工况下总体动能还不够充沛，所以考虑提高进水速度。

图3 进水流速为0.3 m·s-1时的三维和x = 0截面动能分布图

2.1.2 进水速度为0.6 m·s-1的模拟结果

图4是进水流速为0.6 m·s-1的模拟结果，通过该图可见，流速分布依然比较均匀，尽管低位区域的动能分布改善不大，但总体动能更加充沛，加强了整体流场的流动性。

2.1.3 进水速度为1 m·s-1的模拟结果

继续提高进口流速，图5是进水流速为1 m·s-1的模拟结果，通过该图可见，动能比上一个工况更充沛，但分布的均匀性则明显不及。尤其通过x = 0截面，清晰可见在高度约为z=0.5 m和z=0.9 m的区域之间，有明显分层现象，显示尽管消耗了更多能量来搅拌，却反而不利于沼气发酵料液的均匀混合。

图4 进水流速为0.6 m·s-1时的三维和x=0截面动能分布图

图5 进水流速为1 m·s-1时的三维和x = 0截面动能分布图

2.1.4 进水速度为1.5 m·s-1的模拟结果

继续提高进水速度至1.5 m·s-1，通过模拟结果图6可见，速度提高后整个流体区域的动能都非常充沛，但中部的分层现象更加严重。可见，当进水速度达到一定程度后，随着速度的持续增大，流态越来越差。

图6 进水流速为1.5 m·s-1时的三维和x = 0截面动能分布图

2.2 发酵后期的液-液两相流模拟

图7是秸秆原料持续发酵30 d的形态变化，初进入发酵罐的秸秆原料呈较完整的颗粒或长条状，与水相对分离，但经过发酵，会逐渐成为胶融状的假塑性流体[16]。

图7 秸秆原料持续发酵后的形态变化

发酵后期的流态已迥异于纯水，可将秸秆胶融拟流体视为一种粘性流体，则此时的流场为水和秸秆的液-液两相流。可通过模拟计算各流速下，次要相(秸秆)在主相(水)中所占的体积分数在流场中的分布情况，来说明其混合均匀性。

根据实验测定，经过30天充分发酵后，秸秆溶液的总固体(TS)浓度为12.45%，密度为967.56 kg·m-3，运动粘性系数为0.0137 m2·s-1，根据此物性参数设定模拟计算条件。

2.2.1 进水速度为0.6 m·s-1的模拟结果

图8是进水速度为0.6 m·s-1时，第二相(秸秆)所占的体积分数云图。由该图可见，经过水力搅拌达到稳定后，在压力出口(1.1 m高度)的上方区域，形成一个明显的分层，秸秆胶融拟流体所占的体积分数极小，水的体积分数接近100%，大量秸秆原料堆积在下方。此符合沼气发酵后期，清液层浮于料液层上方的经验。这种分层并不利于料液均匀混合，但料液在其下方分布相对均匀，只有进水速度方向上有一个直径不大的通路。可见将压力出口设计得尽量高，可尽量减少清液层所占体积，提高有效容积。

图8 进水流速为0.6 m·s-1时的第二相(秸秆)体积分数三维和x = 0截面分布云图

2.2.2 进水速度为1 m·s-1的模拟结果

图9是进水速度为1 m·s-1的模拟结果。由该图可见，大部分区域秸秆的体积分数分布仍然比较均匀，图5中的分层现象在此时却并未出现。相反，清液层由于受到更大进水速度的冲击而被压缩得更小。可见，在粘性剧增的情况下，适宜的进水速度大大提高了。所以考虑继续提高进水速度，以图将清液层的体积压缩得更小。

2.2.3 进水速度为1.5 m·s-1的模拟结果

图10是进水速度为1.5 m·s-1的模拟结果。秸秆的体积分数分布并不特别均匀，高位的清液层虽然继续变小，但并未被完全打破，相反，在z=0.25 m左右的高度，又形成了一个新的分层，将整个流体区域分成两个浓度差异较大的区间。这可能是多余的动能在中下部形成流体涡旋造成的。可见，过大的进水速度，非但不利于搅拌，反而会形成新的分层，极大降低混合的均匀性，浪费有效容积，不利于提高产气率。

图9 进水流速为1 m·s-1时的第二相(秸秆)体积分数三维和x = 0截面分布云图

图10 进水流速为1.5 m·s-1时的第二相(秸秆)体积分数三维和x = 0截面分布云图

2.3 对比结果讨论

通过不同工况下改变进水流速的模拟计算，沼气料液的搅拌流场得到可视化分析，容易发现无论是清水介质，还是混合两相介质，搅拌中都容易出现明显的分层现象，这是优化搅拌流场首先需要避免的问题。通过高位分散式出口设计和提高底部搅拌的进水流速可以尽量减少清液层所占的体积，但在实际模拟中也发现，太大的进水速度会在罐体中部或底部形成新的分层，不利于均匀混合。可见，在优化沼气工程设计搅拌方案的实际工作中，应该首先使用数值模拟方法，可视化分析搅拌流场形态，确定合适的搅拌参数，尤其是进水流速。

而通过多相流模拟，我们又可以发现，在多相流工况下，适宜的进水流速又发生了改变。以本文所用罐体结构为例，清水介质采用0.6 m·s-1为进水速度，其基本流场优于1 m·s-1，但发酵进行到后期，由于粘性提高，适宜的搅拌速度亦随之提高，1 m·s-1亦不会造成明显分层，则比0.6 m·s-1更加合适，但进水流速继续提高到1.5 m·s-1时，便发现多余的动能会在罐体中下部形成新的明显分层，甚至形成“死区”[17]，对发酵效率非常不利，表明该参数设计偏大，需要优化。可见，在特定的罐体结构中，适宜的进水速度随着粘性增大而增大，并可以通过数值模拟方法定量计算得出。除了罐体结构，死区的形成和原料自身的流变特性也具有密切关系[18]，所以需要在模拟中，根据原料的实际物性特征，正确设置模拟相态参数，获得准确的模拟结果，以供可视化分析。

3 结论

沼气发酵不能完全避免清液层和发酵层的明显分层现象，但可以通过搅拌参数的精心设计尽量避免。进水流速是搅拌设计的关键参数，可以利用CFD方法模拟不同流速的流场优化效果，既避免流速过低不能打破分层，又避免流速过高造成第二相堆积。在本文示例的罐体结构中，适宜的进水速度随着粘性增大而增大，并可以通过CFD数值模拟方法定量计算得出，是现代沼气工程设计的先进方法。