CN107561938B - 微界面强化反应器反应速率构效调控模型建模方法 - Google Patents

Abstract

本发明公开了一种微界面强化反应器反应速率构效调控模型建模方法，基于Levenspiel的理论构建了适用于微界面强化反应器的反应速率构效模型。采用本发明的方法构建的反应速率构效调控模型可以很直观的看出气泡直径、气液传质系数及传质阻力等对反应速率的影响，也把反应体系的气泡直径与反应效率(能效和物效)与体系理化特性、微界面特性、传质特性和反应器结构用数学方法关联起来，从而实现可通过调整结构参数和操作参数以获得反应过程能效物效的最大化目标，或者在给定反应目标(任务)和能耗物耗下，设计出高效的反应器结构。

Description

微界面强化反应器反应速率构效调控模型建模方法

技术领域

本发明属于化工制造、反应器、建模技术领域，具体涉及一种微界面强化反应器反应速率构效调控模型建模方法。

背景技术

氧化、加氢、氯化等多相反应在化工生产过程中广泛存在，其反应速率一般受制于传质过程。气液反应的传质速率主要受液侧(或气侧)传质系数及气液相界面积a共同影响。已有研究表明，a对体积传质系数的影响程度更大，且容易调控。因此，增大a被视为提高受传质控制的气液反应体系反应效率的特别有效的途径。

气泡Sauter平均直径d₃₂是决定a大小的关键参数之一，它们主要受气泡间及气液两相间相互作用力影响。气泡聚并和分裂则分别是上述两种作用力的结果，并影响气泡直径的大小。因此，气泡聚并和破裂作为气泡的介观尺度行为，是决定a大小的深层次原因。关于气泡聚并和分裂行为的研究由来已久，普遍认为能量耗散率和d₃₂是重要的影响因素。事实上，d₃₂能够影响a及体积传质系数大小，是决定气液反应速率的核心因素^[1]。研究显示，当d₃₂逐渐减小时，体积传质速率逐渐增大；特别是当d₃₂小于1mm时，体积传质速率随d₃₂的减小以类似于指数形式较快增大。因此，尽可能地减小d₃₂能够强化气液传质并最终增大反应速率。

鼓泡反应器和搅拌-鼓泡反应器是工业上最传统和常用的气液反应器。如PX氧化制TA的塔式鼓泡反应器，气泡直径通常大于10mm，乃至几厘米级，其传质界面面积十分有限，因此必须将反应器做得很大，以提高反应速率，同时必须通过增加鼓气量来促进液体湍流，使气含率提高，进而增大界面面积，但此举必然降低空气中氧的利用率，增大压缩机功率和尾气排放，导致能耗过渡和物料损失及环境污染。从湍流动力学角度看，传统上用得最广泛的搅拌-鼓泡式气液反应器内大多形成对气泡宏观运动有影响但对气泡破碎作用甚微的大涡，气泡不能有效破碎，故气泡直径偏大，传质面积受限，以致反应效率偏低。为强化气液传质，塔式鼓泡反应器一般增设气体分布板、静态混合器等内件以加强混合，而搅拌釜则需安装不同结构的搅拌桨或内筒等结构，以增加液层的含气量。尽管如此，这两种反应器内的气泡直径通常为5～20mm，所提供的单位体积中的相界面积均十分有限，一般小于100m²/m³,故反应效率不可能获得突破性提高。因此，工业上经常通过高温高压和加大气量来提高气含率和相界面积，但这对反应过程的能耗、物耗及反应选择性都有重大的负面影响。

由于研发气泡的微破碎技术十分重要，故最近10年来，英、美、德、日等国的大学与研究机构开始关注和研发超细气泡技术^[2-11]，但其研究成果具有下列共同缺陷：

(1)采用机械破碎、流体撞击、超声等手段虽可得到一定量微米级尺度的气泡，但气液比(气体体积与液体体积之比)太低，一般低于1％，上限不超过5％。此外，产生微气泡的设备能耗和制造成本也太高。

(2)国内外尚没有基于液相为连续相且高度湍流的微气泡体系特性，提出过系统化的微界面传质强化理论、微气泡测试与表征方法、微界面强化反应器构效调控理论及相关数学模型。

基于上述原因，虽有零星的应用实验结果发表，但尚未有规模化的工业应用报道，尤其在化工制造领域的应用基本上还处于空白状态。

当今时代的化工生产，基于创新、绿色、环保总体考虑，其生存和发展依赖于对材料与过程技术大幅创新。提高反应与分离过程的原子经济性对降低能耗、物耗、增强竞争力至关重要。基于此，我们提出“微界面传质强化反应-精细分离集成系统”新技术，试图从最基本的超细气液颗粒特性研究出发，解决在高度湍流状态下超细颗粒体系中，微界面化学反应器内流体流动、传质、反应、能量转换等构效调控全过程所涉及的理论、技术与应用问题。

本发明涉及的超细气液颗粒指的是超细气泡(或超细液滴)，是颗粒当量直径处在1μm≤d₀＜1mm的微米级气液颗粒。在反应体系中，超细气液颗粒形成了超细界面(或微界面)，超细界面的形成大幅强化了传质与反应速率，尤其是受传质控制的反应体系。

需要强调，经典的气液混合理论一般基于毫米－厘米级气液颗粒特征，目前最为合理的方法是多尺度能量最小化原理(EMMS)^[12]。当前大多研究工作均是针对传统气液反应器内的毫米级上气泡^[13,14]，很少涉及超细颗粒体系。针对超细颗粒体系的混合、传质与反应特性，必须建立新的计算模型、测试与表征方法，以及构效调控模型，为此必须研究新的设备结构、能量输入方式与转换模式，从而形成全新的适合于超细颗粒反应体系的计算软件与硬件平台，为我国的过程工业生产技术上一个新台阶提供技术与装备支持。

气泡流体力学行为是决定气液传质效率的决定性因素，其特征参数一般包括：气泡Sauter平均直径d₃₂、气含率φ_G、气泡上升速度v及滑移速度v_s。对诸如加氢、氧化反应等由液膜扩散控制的气液反应体系而言，这些参数对体系气液相界面积a和传质系数(k_G、k_L)产生重要影响，并最终决定反应速率。为实现气液反应的构-效调控，须从理论上尽可能准确地建立气泡流体力学行为数学模型。

气液多相反应的反应速率不仅取决于本征反应动力学，还极大地受制于体系的气液传质情况^[15]。当反应体系的工艺条件(如催化剂、温度、压力、物料配比等)确定后，该体系的反应速率主要受传质影响。然而，许多研究者往往偏重于对本征反应动力学影响因素的研究而忽略了反应过程传质和传热的研究，其结果往往会导致反应器规模化时出现低效甚至失败。

以过程工业中常见的氧化、加氢反应为例，它们是典型的受传质影响的多相反应。其特征为：(1)反应工艺条件确定后，其本征反应动力学受其它因素影响较小^[16]；(2)反应过程的阻力主要为气液、液固界面的传质阻力^[15]；(3)反应器结构参数、流体运动方式决定体系气-液及液-固界面的传质性能，并最终对反应速率产生决定性影响。

所谓构效调控数学模型，是指将超细气液颗粒反应体系之反应效率(能效和物效)与体系理化特性、微界面特性、传质特性和反应器结构用数学方法关联起来，从而实现可通过调整结构参数和操作参数以获得反应过程能效物效的最大化目标，或者在给定反应目标(任务)和能耗物耗下，设计出高效的反应器结构。而对于微界面强化反应器而言，此方面工作在国际上仍为空白。

参考文献

[1]Levenspiel O.Chemical Reaction Engineering[M].Wiley New York etc.,1972.

[2]Xu JH,Li SW,Chen GG,LuoG..Formation of monodispersemicrobubbles ina microfluidic device[J].AIChE Journal,2006,52(6):2254-2259.

[3]Li P and Tsuge H.Ozone transfer in a new gas-induced contactorwith microbubbles[J].Journal of Chemical Engineering of Japan,2006,39(11):1213-1220.

[4]Muroyama K,Imai K,Oka Y,Hayashi J,Mass transfer properties in abubble column associated with micro-bubble dispersions[J].ChemicalEngineering Science,201,100:464-473.

[5]Maeda Y,Hosokawa S,Baba Y,Tomiyama Akio.Generation mechanism ofmicro-bubbles in a pressurized dissolution method[J].Experimental Thermal andFluid Science,2015,60:201-207.

[6]Hasegawa H,Nagasaka Y,Kataoka H.Electrical potential ofmicrobubble generated by shear flow in pipe with slits.Fluid DynamicsResearch,2008,40(7-8):554-564.

[7]Weber J and Agblevor F.Microbubble fermentation ofTrichodermareesei for cellulase production[J].Process Biochemistry,2005,40(2):669-676.

[8]Rehman F,Medley GJ,Bandulasena H,Zimmerman WB.Fluidic oscillator-mediated microbubble generation to provide cost effective mass transfer andmixing efficiency to the wastewater treatment plants[J].Environmentalresearch,2015,137:32-39.

[9]Stride E and Edirisinghe M.Novel microbubble preparationtechnologies[J].Soft Matter,2008,4(12):2350.

[10]Druzinec D,Salzig De,Kraume M,Czermak P.Micro-bubble aeration inturbulent stirred bioreactors:Coalescence behavior in Pluronic F68 containingcell culture media[J].Chemical Engineering Science,2015,126:160-168.

[11]李宝璋,尚龙安,姜信真.倒喷射式环流反应器的研究[J].西北大学学报(自然科学版).1989,04:65-69.

[12]Chen JH,Yang N,Ge W,Li JH.Stability-driven structure evolution:exploring the intrinsic similarity between gas-solid and gas-liquid systems[J].Chinese Journal of Chemical Engineering.2012,20(1):167-177.

[13]Hinze JO.Fundamentals of the hydrodynamic mechanism of splittingin dispersion processes[J].AIChE Journal.1955,1(3):289-295.

[14]Zhong S,Zou X,Zhang ZB,Tian HZ.A flexible image analysis methodfor measuring bubble parameters[J].Chemical Engineering Science,2016,141(17):143-153.

[15]Levenspiel O.Chemical reaction engineering[M].Wiley New Yorketc.,1972.

[16]Deckwer W,Schumpe A.Improved tools for bubble column reactordesign and scale-up[J].Chemical Engineering Science.1993,48(5):889-911.

发明内容

本发明的目的在于，克服现有技术的缺陷，提供一种微界面强化反应器宏观速率构效调控模型建模方法。本发明基于对Levenspiel所认为的多项体系的宏观速率表达式的研究，通过合理的假设和简化，建立一个应用于微界面反应体系的宏观速率表达式。并且通过对各个参数的定义和推导，形成一套完整的理论体系，来构建宏观速率构效调控模型，以此来为指导微界面强化反应器的设计和调控。

为实现上述技术目的，本发明采用如下技术方案：

获取微界面强化反应器反应体系工艺条件，包括催化剂、温度、压力、物料配比；计算气液反应反应速率；基于Levenspiel的理论，多相体系的反应速率由下式表示：

其中，r_A为反应物A的反应速率；H_A为反应物A的亨利常数；k_G、k_L、k_S分别为气侧、液侧及液固传质系数；a、a_S分别为气液相界面积、液固相界面积；k_A为基于反应速率r_A的一级本征反应速率常数；

为催化剂颗粒内部溶剂B的平均浓度；χ_A为反应速率常数k_A时的一级本征反应有效因子，表征因催化剂孔扩散导致的反应速率降低的程度；f_s为催化剂装载率；P_G为气泡内压；

选定χ_A＝1；

是催化剂表面溶剂的浓度，在数量上与f_s形成反向对应关系，即有：

忽略式(1)分母第三项，将气液反应反应速率方程简化为：

P_G表示为：

式中，ρ_L为液体密度；H₀为反应器内初始液位高度；d₃₂为反应器内气泡Sauter平均直径；v₃₂为与d₃₂对应的反应器内气泡平均上升速度；v_G为反应器内的表观气速；σ_L为液体表面张力；

将(4)带入(3)，得到式(5)：

进一步的，所述传质系数采用如下方式获取：

气侧传质系数k_G：根据膜理论，气侧传质系数k_G大小由气相扩散系数D_G及有效气膜厚度δ_G决定；

基于Chapman-Enskog动理论，D_G的一般形式为：

式中，T为气体温度；M_A、M_B分别为气体A及溶剂B的摩尔质量；P_G为反应器内气泡内气体平均压力；采用Gedde静止球模型建立k_G理论计算模型：

式中，d₃₂为反应器内气泡Sauter平均直径；t₃₂为直径为d₃₂的气泡在体系中的停留时间，其表达式为：

式中，H₀为反应器内初始液位高度；v₃₂为与d₃₂对应的反应器内气泡平均上升速度；v_G为反应器内的表观气速，Q_G为反应器内通气体积流量；

将式(8)代入上式(7)后得：

液侧传质系数k_L：建立基于表面更新理论的渗透模型，具体如下：

式中t_c为液体微元暴露时间；D_L为液相扩散系数；

依据速度滑移模型的定义，将t_c利用式(11)定义：

t_c＝d₃₂/v_S (11)

式中v_s为气泡滑移速度；

将式(11)代入式(10)得：

根据气泡滑移速度v_s的定义：

v_S＝v₃₂-u (13)

式中u为气泡间流体速度；

当液体循环流量为Q_L，对于均匀上升气流，u的表达式为：

式中，v_L为表观液速，

S₀为反应器横截面积，D₀为反应器直径；φ_G为反应器内气含率，假定为φ_G＝v_G/v₃₂；

将式(14)代入式(13)得：

将式(15)代入式(12)得：

式中，液相扩散系数D_L基于修正的Wilke–Chang方程获取，如式(17)所示：

进一步的，所述气液相界面积a采用如下方式获取：

步骤1：微界面强化反应器气液相界面积通用表达式建立：

假设微界面强化反应器内气泡为刚性球体，设定气泡最大和最小直径分别为d_max和d_min；当反应器稳态工作时，有：

式中，a为反应器内气液相界面积，即单位体积气液混合物中气泡的表面积；x、N、F(x)分别为气泡直径、稳态时反应器内气泡总数及气泡粒径分布概率密度函数；V_T为气液混合物总体积；

设反应器内初始液位高度和反应器直径分别为H₀和D₀，体系气含率为φ_G，计算V_T和反应器横截面积S₀：

将式(19)代入式(18)，得：

依据φ_G的定义：

将式(19)代入式(22)，得：

在稳态情况下，将慢反应体系中气泡单程上升时由于反应消耗的体积及数量忽略不计，则进出反应器体系的气体分布基本相等，且气泡逃逸概率f(x)为：

f(x)NF(x)＝n_II(x) (24)

上式中，n_I和I(x)分别为单位时间进入反应器的气泡数量和粒径分布概率密度函数；当反应体系的通气量为Q_G，依据气体体积守恒原则，有：

将式(25)代入(24)并积分，得：

将式(26)代入式(23)后，得：

式(27)中，v_G为表观气速，由下式表示：

v_G＝Q_G/S₀ (28)

假设气泡逸出概率f(x)等于气泡在体系中停留时间的倒数，则有：

f(x)＝v(x)(1-φ_G)/H₀ (29)

式中v(x)是直径为x的气泡上升速度；将式(29)分别代入式(26)、(27)，得：

将式(30)代入式(21)后，得：

假设气泡进入反应器后一般不发生聚并和破裂现象；且由于是慢反应过程，气泡从反应器底部上升至顶部的单程运动过程中，所反应消耗的量可忽略不计，即稳态时反应器内各处的气泡粒径分布与初始气泡分布一致，也就是：

F(x)＝I(x) (32)

将式(32)代入式(31)并进行离散化处理后得到a的计算公式为：

步骤2：微界面强化反应器气含率通用表达式建立：

基于已有a的另一通用表达式，即式(34)：

a＝6φ_G/d₃₂ (34)

式中，d₃₂为气泡Sauter平均直径；对比(33)和(34)，当φ_G定义为下式(35)时，上述两式等价：

式(35)即为φ_G通用计算公式；

步骤3：获取微界面强化反应器气液体系内气泡上升速度v₃₂；

φ_G是体系宏观特征参数，与d₃₂类似，假定：

φ_G＝v_G/v₃₂ (36)

可得：

对于微米级气泡，考虑升力、气泡周围压力梯度、气泡对其周围湍流场响应时间τ_b(s)对曳力的影响，微气泡运动的矢量方程可表示为式(38)：

式中，u为气泡间流体速度；y(t)表示t时刻气泡的位置；g为重力加速度；ω为气泡涡量；τ_b为气泡响应时间；

对于微界面强化反应器的微气泡体系，反应器内为无旋势流场，因此ω＝0；假设气泡进入体系后宏观上为匀速上升，并且稳态时，反应器内湍流各向同性，气泡稳定上升，则Du/Dt＝0；dv₃₂/dt＝0；气泡响应时间τ_b表示为：

式中，v₀为无限大静止液相中单个气泡上升速度；

仅考虑反应器内气泡一维运动情况并结合以上假设，将式(38)简化：

v₃₂＝v₀±u (40)

分析反应器内气泡周围液体的三种运动情况：

a.气泡周围液体宏观竖直向下运动；

此时：

v₃₂＝v₀-u (41)

当体系中气泡较小时，基于反应器工业生产需求，须满足：

v₀≥u (42)

结合式(14)(41)，有：

且：

对式(43)化简后得：

v₃₂ ²+(v_L-v_G-v₀)v₃₂+v_Gv₀＝0 (45)

求解式(45)，获得唯一实解：

b.气泡周围液体宏观运动速度为零；

根据微界面强化反应器中微气泡主导的反应体系，将Fan通过实验拟合获取的v₀的表达式化简得：

式中，σ_L为液体表面张力；μ_L为液体动力粘度；ρ_L为液体密度；

c.气泡周围液体宏观竖直向上运动：

此时：

v₃₂＝v₀+u (48)

并且：

v₃₂ ²-(v_L+v_G+v₀)v₃₂+v_Gv₀＝0 (49)

整理后得：

根据所述a、b、c，气泡上升速度v₃₂依据液体宏观流动方向表示为：

式中，v_L前的“+”表示气泡周围液体宏观竖直向上运动的情况，“-”则表示气泡周围液体宏观竖直向下运动的情况；上式中v₀依据式(42)计算。

进一步的，所述反应器内气泡Sauter平均直径d₃₂采用如下方式获取：

(1)以微界面强化反应器最大气泡直径d_max和最小气泡直径d_min为自变量，气泡Sauter平均直径d₃₂为因变量，建立d_max、d_min和d₃₂间的关系；具体步骤如下：

设x,m,n分别为反应器气液体系中的气泡粒径、气泡粒径几何自然对数的均值和标准差，获取气泡粒径x的概率密度函数：

气泡粒径满足此分布时的气泡Sauter平均直径d₃₂为：

d₃₂＝exp(m+2.5n²) (53)

气泡粒径x呈对数正态分布，因此lnx的数学期望(算术平均值)为：

根据气泡粒径x的概率密度函数绘制气泡粒径概率密度图，当

时，概率密度最大；即此处的一阶导数为0：

将方程(54)代入(54)得到方程(55)：

由(54)、(55)可得：

由于：

将方程(52)代入(57)并化简后可得：

令：

则上式简化为：

方程(59)左端为误差函数，与标准误差函数的差别在于积分限的不同，将式(56)分别代入上述积分上下限，并将方程(59)转化为标准误差函数后可得：

方程(60)中，erf(·)为误差函数；

对于如下形式的误差函数：

采用泰勒级数展开进行近似计算，泰勒级数展开依据误差函数自变量的取值范围不同而采用不同的形式，当z≤4时，erf(z)可展开为：

由于：

当d_max/d_min为1000时：

而根据方程(62)：

因此，当：

即：

时，方程(60)近似成立；

此外，方程(60)成立的条件与n及d_max/d_min的大小有关，且n受d_max/d_min的大小的制约；构建气泡粒径累积概率密度g(n)以考察n和d_max/d_min对方程(60)成立条件的影响，令气泡粒径累积概率密度g(n)为：

绘制g(n)～n关系曲线；获取确保方程(60)成立的n的可取值范围与d_max/d_min的关系；

取不等式(67)的等号条件，即：

由式(56)及(69)确定m和n，进而由方程(53)建立d₃₂基本数学模型；其结果如下：

(2)基于Kolmogorov-Hinze理论，构建微界面强化反应器最大气泡直径d_max、最小气泡直径d_min与反应器参数间的关系；

能使气泡破裂的最小湍流涡尺度是Kolmogorov尺度的11.4～31.4倍，假设此倍率为11.4，由于湍流涡仅能破碎直径大于其尺度的气泡，因此，气泡直径最小值d_min与该湍流涡尺度一致，即：

d_min＝11.4(μ_L/ρ_L)^0.75ε^-0.25 (72)

基于Kolmogorov-Hinze理论，最大气泡直径d_max由下式(73)确定：

d_max＝ε^-2/5(σ_LWe_crit/2ρ_L)^3/5 (73)

其中，ε为能量耗散率；μ_L为液体动力粘度；σ_L为液体表面张力；ρ_L为液体密度；We_crit为气泡破碎临界weber数；

基于气泡破碎的共振理论确定We_crit：

其中，α₂为气泡体积模量，α₂＝2,3,...；当α₂越大，气泡高阶振动越激烈，气泡就越小，对于超细气泡颗粒选择α₂＝2，即We_crit＝1.24；

此时：

d_max＝0.75(σ_L/ρ_L)^0.6ε^-0.4 (75)

将式(72)、式(75)代入式(71)求解d₃₂。

进一步的，所述能量耗散率ε采用如下方式获取：

步骤100：将微界面强化反应器总的能量耗散率ε的计算划分为微界面强化反应器内三个不同区域能量耗散率的总和，包括反应器主体区鼓泡区的能量耗散率ε_R，气液破碎区的ε_mix以及气液出口区的ε_pl；

步骤110：其中，反应器主体区鼓泡区的能量耗散率ε_R采用如下方式计算：

反应器气体鼓泡过程中，依据气泡对体系做功，ε_R表示为：

其中，Q_G为反应器内通气体积流量，m³/s；S₀为反应器横截面积，m²；

步骤120：计算气液破碎区的ε_mix：

基于ε_mix经典计算模型，假设气液混合为绝热过程并忽略液体势能变化，忽略气体质量流量，并使能量耗散率的单位统一为W/Kg，获取计算公式如下：

其中，L_mix为气泡破碎区长度，m；P₀、P₁分别为气泡破碎区入口液体静压及出口气液混合物压力，Pa；λ₁为气液体积流量之比：K₁为喷嘴直径与气泡破碎区直径的比值，K₁＝D_N/D₁；S₁为气泡破碎器横截面积，m²；ρ_L为液体密度，kg/m³；Q_L为反应器内液体循环体积流量，m³/s；

λ₁＝Q_G/Q_L (78)

步骤121：计算气泡破碎区入口液体静压P₀及出口气液混合物压力P₁：

忽略气泡破碎区管壁摩擦损失，则：

其中，φ_mix为气泡破碎区气含率，按下式计算：

忽略管道摩擦及喷嘴处能量损失，依据能量守恒原理，体系实际获得的能量E₀为：

即：

由式(79)、(82)得：

步骤122：计算气泡破碎区长度L_mix：

L_mix通过测量气液破碎区管内壁压力突变来确定，或通过如下方式确定：

其中：P_H为气液破碎区上方气压，Pa；ρ_MZ为气液破碎区内气液混合物密度，Kg/m³；v_N为射流口的射流速度，m/s；U_e,max为气液破碎区涡旋最大返回速度，m/s；

P_H由伯努利方程推得：

P_H≈P_G0 (85)

式中，P_G0为供气压力，Pa；

ρ_MZ通过下式计算：

ρ_MZ＝ρ_Gφ_mix+ρ_L(1-φ_mix)≈ρ_L(1-φ_mix) (86)

式中，ρ_G为气体密度，g/m³；

考虑气液破碎区气液混合物流速的影响，U_e,max为射流口射流流速与气液破碎区气液混合物流速的矢量合成的结果，采用下式计算U_e,max：

将式(85)(87)代入式(84)，并化简后可得：

获取反应器气泡破碎器长度L_b，并根据式(88)计算L_mix；

①当L_mix＜L_b时，以式(88)的计算结果为L_mix的实际数值；

②当L_mix≥L_b时，说明射流能量近似完全在气泡破碎器区域耗，则：

L_mix＝L_b (89)

步骤130：计算气液出口区的ε_pl；

假设气泡在气液出口区呈均匀分布状态，此区域的能量耗散速率ε_pl由下式计算：

反应器结构设计时保证λ₁可调范围足够大，通过实验确定反应器基本结构参数间的关系为K₁＝0.5，L_b＝13D₁；代入前述相应表达式并化简可得：

步骤200：确定ε_R、ε_mix和ε_pl各自的数值大小；

步骤210：依据进入反应器的气液流量等于气泡破碎区出口的气液流量平衡原理，得到：

式中，C_L为基于反应器内有效体积πD₀ ²H₀/4的液体循环倍数，即每小时液体循环总体积与反应器有效体积的比值；u₁为气泡破碎器出口气液混合物线速度，m/s；λ₁取值0.1～0.5；

由式(94)可知：

则u₁增大时，反应器横截面积S₀也增大；结合式(76)可知，此时ε_R减小；为对反应器不同区域的能量耗散率进行比较，假设：u₁＝3.0m/s；C_L＝20；H₀＝1.5m；由式(94)可得，当λ₁＝0.1～0.5时：

D₀≈19D₁ (95)

选定D₁数值，计算并比较反应器不同区域在不同喷嘴液速下的能量耗散率，确定与气液破碎区的能量耗散率ε_mix相比，反应器主体区鼓泡区、气液出口区的能量耗散率可忽略不计，即ε_mix≈ε；则整个反应器的能量耗散率ε与反应器结构参数之间的数学关系，可由式(77)计算确定，即：

本发明的另一目的在于提供上述方法建立的反应速率构效调控模型，具体如下：

d_min＝11.4(μ_L/ρ_L)^0.75ε^-0.25 (72)

d_max＝0.75(σ_L/ρ_L)^0.6ε^-0.4 (75)

式中，k_G为气测传质系数；k_L为气测传质系数；k_A为基于反应速率r_A的一级本征反应速率常数；D_G为气相扩散系数；D_L为液相扩散系数；d₃₂为反应器内气泡Sauter平均直径；d_min、d_max分别为反应器中最大气泡和最小气泡直径；v₃₂为与d₃₂对应的反应器内气泡平均上升速度；v_G为反应器内的表观气速；v_L为反应器内的表观液速；v₀：无限大静止液相中单个气泡上升速度；S₀：反应器横截面积；D₀：反应器直径；Q_L为反应器内液体循环体积流量；L_mix为气泡破碎区长度；D₁为气泡破碎管直径；λ₁为气液体积流量之比，λ₁＝Q_G/Q_L；Q_G为反应器内通气体积流量；P₀为气泡破碎器入口处液体的静压；P₁为气泡破碎区出口气液混合物压力；ε为能量耗散率；μ_L为液体动力粘度；σ_L为液体表面张力；ρ_L为液体密度。

本发明的又一目的在于提供上述方法在反应器设计中的应用。

本发明的方法适用于微界面强化反应器，其核心在于气泡破碎器。气泡破碎器的原理是高速射流所携带的气体相互撞击进行能量传递，使气泡破碎，其结构参数有L_b、D₁，详细结构见附图1，除此之外该反应器的其它结构参数有D₀、H₀，具体反应器结构相关内容已经公布于发明人在先申请的专利CN106187660A中，本发明中不再赘述。

本发明的方法具有如下有益效果：

(1)采用本发明的方法构建的反应速率构效调控模型可以很直观的看出气泡直径、气液传质系数及传质阻力等对反应速率的影响，也把反应体系的气泡直径与反应效率(能效和物效)与体系理化特性、微界面特性、传质特性和反应器结构用数学方法关联起来，从而实现可通过调整结构参数和操作参数以获得反应过程能效物效的最大化目标，或者在给定反应目标(任务)和能耗物耗下，设计出高效的反应器结构。而对于微界面强化反应器而言，此方面工作在国际上仍为空白。

(2)本发明的方法构建了d_max、d_min和d₃₂的直接计算关系，而不再采用实验拟合的方式获取d₃₂的具体数值，极大减少了在反应器中应用时产生的误差；现有的方法构建的d₃₂模型多针对的是鼓泡反应器(Bubble column,BC)和鼓泡搅拌釜反应器(Bubbling-stirredreactor,BSR)内的气液体系，或射流泵(gas-liquid jet bump,GLJB)内的空气-水体系。而对于工业微界面反应器(MIR)，则不一定适用，其原因在于：①MIR内气泡破碎机制有别于上述反应器；②工业气液反应体系中可能涉及到高粘度液体，而现有技术并没有考虑液体粘度对d₃₂的影响；而本发明的方法适用于工业微界面反应器(MIR)，其通用性更好。

(3)本发明的方法还构建了能量耗散率ε的表达式，ε是决定气泡破裂机制的关键参数。而ε的增大，意味着气泡更容易破裂，即可以得到尺度更小的气泡，但是针对微米级尺度的能量耗散率的表达式，目前相关研究还是空白，本发明的方法通过严谨的推导过程实现了微界面强化反应器能量耗散率ε的数学模型的构建，为指导新型反应器的设计奠定理论基础。

附图说明

图1是一种微界面强化反应器结构示意图，用于说明本发明建模方法在反应器装置中的应用；其中1-反应器，2-泵前阀，3-循环泵，4-泵后阀，5-液体流量计，6-换热器，7-气泡破碎器，8-测温仪，9-下降管，10-气体阀，11-气体流量计，12-气相入口，13-压力表，14-液位计；D₀-反应器直径，H₀-反应器内初始液位高度，D₁-气泡破碎管直径，L_b-气泡破碎区长度。

具体实施方式

实施例1

本实施例具体说明本发明建模方法的技术方案。

本发明的方法具体包括如下步骤：

获取微界面强化反应器反应体系工艺条件，包括催化剂、温度、压力、物料配比；计算气液反应反应速率；

基于Levenspiel的理论，多相体系的反应速率由下式表示：

上式中，r_A为反应物A的反应速率，mol(A)/m^-3(reactor)·s；H_A为反应物A的亨利常数，Pa·m³/mol；k_G为气侧传质系数，mol/(Pa·m²·s)、k_L、k_S分别为液侧及液固传质系数，m/s。a、a_S分别为气液相界面积、液固相界面积，m²/m³；k_A为基于反应速率r_A的一级本征反应速率常数，m³(A)/m³(cat)·s；

为催化剂颗粒内部溶剂B的平均浓度，mol/l；χ_A为反应速率常数k_A时的一级本征反应有效因子，表征因催化剂孔扩散导致的反应速率降低的程度；f_s为催化剂装载率，m³(cat)/m³(reactor)；P_G为气泡内压，Pa。

基于反应动力学，可以认为上式中分母各项均为反应阻力。随着气液反应体系高性能催化剂的出现(一般以分子态或溶液态存在于反应体系中)，可以认为气液反应多在拟均相催化剂下发生。而且由于催化剂的装载率f_s相对反应器体积非常小(一般为1/1000～1/10000),故可认为气液催化反应发生于催化剂外表面，故χ_A＝1。又由于

是催化剂表面溶剂的浓度,其数值很大，它在数量上与f_s形成反向对应关系，即有：

此外，相对于气液传质速率而言，k_Sa_S要大很多。因此，上式分母第三项可以忽略。因此气液反应反应速率方程可简化为：

上式反映出传质和本征反应速率常数对反应速率的影响。

其中，P_G表示为：

式中，ρ_L为液体密度；H₀为反应器内初始液位高度；d₃₂为反应器内气泡Sauter平均直径；v₃₂为与d₃₂对应的反应器内气泡平均上升速度；v_G为反应器内的表观气速；σ_L为液体表面张力；

将(4)带入(3)，得到式(5)：

其中，所述传质系数采用如下方式获取：

气侧传质系数k_G：

根据膜理论，气侧传质系数k_G大小由气相扩散系数D_G及有效气膜厚度δ_G决定；

基于Chapman-Enskog动理论，D_G的一般形式为：

式中，T为气体温度；M_A、M_B分别为气体A及溶剂B的摩尔质量；P_G为反应器内气泡内气体平均压力；采用Gedde静止球模型建立k_G理论计算模型：

式中，d₃₂为反应器内气泡Sauter平均直径；t₃₂为直径为d₃₂的气泡在体系中的停留时间，其表达式为：

式中，H₀为反应器内初始液位高度；v₃₂为与d₃₂对应的反应器内气泡平均上升速度；v_G为反应器内的表观气速，Q_G为反应器内通气体积流量；

将式(8)代入上式(7)后得：

液侧传质系数k_L：

建立基于表面更新理论的渗透模型，具体如下：

式中t_c为液体微元暴露时间；D_L为液相扩散系数；

依据速度滑移模型的定义，将t_c利用式(11)定义：

t_c＝d₃₂/v_S (11)

式中v_s为气泡滑移速度；

将式(11)代入式(10)得：

根据气泡滑移速度v_s的定义：

v_S＝v₃₂-u (13)

式中u为气泡间流体速度；

当液体循环流量为Q_L，对于均匀上升气流，u的表达式为：

式中，v_L为表观液速，

S₀为反应器横截面积，D₀为反应器直径；φ_G为反应器内气含率，假定为φ_G＝v_G/v₃₂；

将式(14)代入式(13)得：

将式(15)代入式(12)得：

式中，液相扩散系数D_L基于修正的Wilke–Chang方程获取，如式(17)所示：

所述气液相界面积a采用如下方式获取：

步骤1：微界面强化反应器气液相界面积通用表达式建立：

假设微界面强化反应器内气泡为刚性球体，设定气泡最大和最小直径分别为d_max和d_min；当反应器稳态工作时，有：

式中，a为反应器内气液相界面积，即单位体积气液混合物中气泡的表面积；x、N、F(x)分别为气泡直径、稳态时反应器内气泡总数及气泡粒径分布概率密度函数；V_T为气液混合物总体积；

设反应器内初始液位高度和反应器直径分别为H₀和D₀，体系气含率为φ_G，计算V_T和反应器横截面积S₀：

将式(19)代入式(18)，得：

依据φ_G的定义：

将式(19)代入式(22)，得：

在稳态情况下，将慢反应体系中气泡单程上升时由于反应消耗的体积及数量忽略不计，则进出反应器体系的气体分布基本相等，且气泡逃逸概率f(x)为：

上式中，n_I和I(x)分别为单位时间进入反应器的气泡数量和粒径分布概率密度函数；当反应体系的通气量为Q_G，依据气体体积守恒原则，有：

将式(25)代入(24)并积分，得：

将式(26)代入式(23)后，得：

式(27)中，v_G为表观气速，由下式表示：

v_G＝Q_G/S₀ (28)

假设气泡逸出概率f(x)等于气泡在体系中停留时间的倒数，则有：

f(x)＝v(x)(1-φ_G)/H₀ (29)

式中v(x)是直径为x的气泡上升速度；将式(29)分别代入式(26)、(27)，得：

将式(30)代入式(21)后，得：

假设气泡进入反应器后一般不发生聚并和破裂现象；且由于是慢反应过程，气泡从反应器底部上升至顶部的单程运动过程中，所反应消耗的量可忽略不计，即稳态时反应器内各处的气泡粒径分布与初始气泡分布一致，也就是：

F(x)＝I(x) (32)

将式(32)代入式(31)并进行离散化处理后得到a的计算公式为：

步骤2：微界面强化反应器气含率通用表达式建立：

基于已有a的另一通用表达式，即式(34)：

a＝6φ_G/d₃₂ (34)

式中，d₃₂为气泡Sauter平均直径；对比(33)和(34)，当φ_G定义为下式(35)时，上述两式等价：

式(35)即为φ_G通用计算公式；

步骤3：获取微界面强化反应器气液体系内气泡上升速度v₃₂；

φ_G是体系宏观特征参数，与d₃₂类似，假定：

φ_G＝v_G/v₃₂ (36)

可得：

对于微米级气泡，考虑升力、气泡周围压力梯度、气泡对其周围湍流场响应时间τ_b(s)对曳力的影响，微气泡运动的矢量方程可表示为式(38)：

式中，u为气泡间流体速度；y(t)表示t时刻气泡的位置；g为重力加速度；ω为气泡涡量；τ_b为气泡响应时间；

对于微界面强化反应器的微气泡体系，反应器内为无旋势流场，因此ω＝0；假设气泡进入体系后宏观上为匀速上升，并且稳态时，反应器内湍流各向同性，气泡稳定上升，则Du/Dt＝0；dv₃₂/dt＝0；气泡响应时间τ_b表示为：

式中，v₀为无限大静止液相中单个气泡上升速度；

仅考虑反应器内气泡一维运动情况并结合以上假设，将式(38)简化：

v₃₂＝v₀±u (40)

分析反应器内气泡周围液体的三种运动情况：

a.气泡周围液体宏观竖直向下运动；

此时：

v₃₂＝v₀-u (41)

当体系中气泡较小时，基于反应器工业生产需求，须满足：

v₀≥u (42)

结合式(14)(41)，有：

且：

对式(43)化简后得：

v₃₂ ²+(v_L-v_G-v₀)v₃₂+v_Gv₀＝0 (45)

求解式(45)，获得唯一实解：

b.气泡周围液体宏观运动速度为零；

根据微界面强化反应器中微气泡主导的反应体系，将Fan通过实验拟合获取的v₀的表达式化简得：

式中，σ_L为液体表面张力；μ_L为液体动力粘度；ρ_L为液体密度；

c.气泡周围液体宏观竖直向上运动：

此时：

v₃₂＝v₀+u (48)

并且：

v₃₂ ²-(v_L+v_G+v₀)v₃₂+v_Gv₀＝0 (49)

整理后得：

根据所述a、b、c，气泡上升速度v₃₂依据液体宏观流动方向表示为：

式中，v_L前的“+”表示气泡周围液体宏观竖直向上运动的情况，“-”则表示气泡周围液体宏观竖直向下运动的情况；上式中v₀依据式(42)计算。

现有d₃₂的具体数值采用实验拟合的方式获取，在反应器中应用时可能产生很大的误差，本发明设计的d₃₂算法构建了d_max、d_min和d₃₂的直接计算关系，极大地减小了在反应器中应用时的误差。

所述反应器内气泡Sauter平均直径d₃₂采用如下方式获取：

(1)以微界面强化反应器最大气泡直径d_max和最小气泡直径d_min为自变量，气泡Sauter平均直径d₃₂为因变量，建立d_max、d_min和d₃₂间的关系；具体步骤如下：

设x,m,n分别为反应器气液体系中的气泡粒径、气泡粒径几何自然对数的均值和标准差，获取气泡粒径x的概率密度函数：

气泡粒径满足此分布时的气泡Sauter平均直径d₃₂为：

d₃₂＝exp(m+2.5n²) (53)

气泡粒径x呈对数正态分布，因此lnx的数学期望(算术平均值)为：

根据气泡粒径x的概率密度函数绘制气泡粒径概率密度图，当

时，概率密度最大；即此处的一阶导数为0：

将方程(54)代入(54)得到方程(55)：

由(54)、(55)可得：

由于：

将方程(52)代入(57)并化简后可得：

令：

则上式简化为：

方程(59)左端为误差函数，与标准误差函数的差别在于积分限的不同，将式(56)分别代入上述积分上下限，并将方程(59)转化为标准误差函数后可得：

方程(60)中，erf(·)为误差函数；

对于如下形式的误差函数：

采用泰勒级数展开进行近似计算，泰勒级数展开依据误差函数自变量的取值范围不同而采用不同的形式，当z≤4时，erf(z)可展开为：

由于：

当d_max/d_min为1000时：

而根据方程(62)：

因此，当：

即：

时，方程(60)近似成立；

此外，方程(60)成立的条件与n及d_max/d_min的大小有关，且n受d_max/d_min的大小的制约；构建气泡粒径累积概率密度g(n)以考察n和d_max/d_min对方程(60)成立条件的影响，令气泡粒径累积概率密度g(n)为：

绘制g(n)～n关系曲线；获取确保方程(60)成立的n的可取值范围与d_max/d_min的关系；

取不等式(67)的等号条件，即：

由式(56)及(69)确定m和n，进而由方程(53)建立d₃₂基本数学模型；其结果如下：

(2)基于Kolmogorov-Hinze理论，构建微界面强化反应器最大气泡直径d_max、最小气泡直径d_min与反应器参数间的关系；

能使气泡破裂的最小湍流涡尺度是Kolmogorov尺度的11.4～31.4倍，假设此倍率为11.4，由于湍流涡仅能破碎直径大于其尺度的气泡，因此，气泡直径最小值d_min与该湍流涡尺度一致，即：

d_min＝11.4(μ_L/ρ_L)^0.75ε^-0.25 (72)

基于Kolmogorov-Hinze理论，最大气泡直径d_max由下式(73)确定：

d_max＝ε^-2/5(σ_LWe_crit/2ρ_L)^3/5 (73)

其中，ε为能量耗散率；μ_L为液体动力粘度；σ_L为液体表面张力；ρ_L为液体密度；We_crit为气泡破碎临界weber数；

基于气泡破碎的共振理论确定We_crit：

其中，α₂为气泡体积模量，α₂＝2,3,...；当α₂越大，气泡高阶振动越激烈，气泡就越小，对于超细气泡颗粒选择α₂＝2，即We_crit＝1.24；

此时：

d_max＝0.75(σ_L/ρ_L)^0.6ε^-0.4 (75)

将式(72)、式(75)代入式(71)求解d₃₂。

其中，所述能量耗散率ε是决定气泡破裂机制的关键参数。本发明基于微界面强化反应器，设计了更适用于微气泡体系的ε计算方法。

所述能量耗散率ε采用如下方式获取：

步骤100：将微界面强化反应器总的能量耗散率ε的计算划分为微界面强化反应器内三个不同区域能量耗散率的总和，包括反应器主体区鼓泡区的能量耗散率ε_R，气液破碎区的ε_mix以及气液出口区的ε_pl；

步骤110：其中，反应器主体区鼓泡区的能量耗散率ε_R采用如下方式计算：

反应器气体鼓泡过程中，依据气泡对体系做功，ε_R表示为：

其中，Q_G为反应器内通气体积流量，m³/s；S₀为反应器横截面积，m²；

步骤120：计算气液破碎区的ε_mix：

基于ε_mix经典计算模型，假设气液混合为绝热过程并忽略液体势能变化，忽略气体质量流量，并使能量耗散率的单位统一为W/Kg，获取计算公式如下：

其中，L_mix为气泡破碎区长度，m；P₀、P₁分别为气泡破碎区入口液体静压及出口气液混合物压力，Pa；λ₁为气液体积流量之比：K₁为喷嘴直径与气泡破碎区直径的比值，K₁＝D_N/D₁；S₁为气泡破碎器横截面积，m²；ρ_L为液体密度，kg/m³；Q_L为反应器内液体循环体积流量，m³/s；

λ₁＝Q_G/Q_L (78)

步骤121：计算气泡破碎区入口液体静压P₀及出口气液混合物压力P₁：

忽略气泡破碎区管壁摩擦损失，则：

其中，φ_mix为气泡破碎区气含率，按下式计算：

忽略管道摩擦及喷嘴处能量损失，依据能量守恒原理，体系实际获得的能量E₀为：

即：

由式(79)、(82)得：

步骤122：计算气泡破碎区长度L_mix：

L_mix通过测量气液破碎区管内壁压力突变来确定，或通过如下方式确定：

其中：P_H为气液破碎区上方气压，Pa；ρ_MZ为气液破碎区内气液混合物密度，Kg/m³；v_N为射流口的射流速度，m/s；U_e,max为气液破碎区涡旋最大返回速度，m/s；

P_H由伯努利方程推得：

P_H≈P_G0 (85)

式中，P_G0为供气压力，Pa；

ρMZ通过下式计算：

ρ_MZ＝ρ_Gφ_mix+ρ_L(1-φ_mix)≈ρ_L(1-φ_mix) (86)

式中，ρ_G为气体密度，g/m³；

考虑气液破碎区气液混合物流速的影响，U_e,max为射流口射流流速与气液破碎区气液混合物流速的矢量合成的结果，采用下式计算U_e,max：

将式(85)(87)代入式(84)，并化简后可得：

获取反应器气泡破碎器长度L_b，并根据式(88)计算L_mix；

①当L_mix＜L_b时，以式(88)的计算结果为L_mix的实际数值；

②当L_mix≥L_b时，说明射流能量近似完全在气泡破碎器区域耗，则：

L_mix＝L_b (89)

步骤130：计算气液出口区的ε_pl；

假设气泡在气液出口区呈均匀分布状态，此区域的能量耗散速率ε_pl由下式计算：

反应器结构设计时保证λ₁可调范围足够大，通过实验确定反应器基本结构参数间的关系为K₁＝0.5，L_b＝13D₁；代入前述相应表达式并化简可得：

步骤200：确定ε_R、ε_mix和ε_pl各自的数值大小；

步骤210：依据进入反应器的气液流量等于气泡破碎区出口的气液流量平衡原理，得到：

式中，C_L为基于反应器内有效体积πD₀ ²H₀/4的液体循环倍数，即每小时液体循环总体积与反应器有效体积的比值；u₁为气泡破碎器出口气液混合物线速度，m/s；λ₁取值0.1～0.5；

由式(94)可知：

则u₁增大时，反应器横截面积S₀也增大；结合式(76)可知，此时ε_R减小；为对反应器不同区域的能量耗散率进行比较，假设：u₁＝3.0m/s；C_L＝20；H₀＝1.5m；由式(94)可得，当λ₁＝0.1～0.5时：

D₀≈19D₁ (95)

选定D₁数值，计算并比较反应器不同区域在不同喷嘴液速下的能量耗散率，确定与气液破碎区的能量耗散率ε_mix相比，反应器主体区鼓泡区、气液出口区的能量耗散率可忽略不计，即ε_mix≈ε；则整个反应器的能量耗散率ε与反应器结构参数之间的数学关系，可由式(77)计算确定，即：

实施例2

本实施例以图1所示的反应器为例，说明实施例1所述建模方法构建的模型在二氧化碳和水体系反应器中的应用。图1的反应器结构可为现有反应器的结构，仅采用本发明的方法进行参数设计，本发明中对反应器的结构不再赘述。

根据实施例1构建的反应速率构效调控模型如下：

d_min＝11.4(μ_L/ρ_L)^0.75ε^-0.25 (72)

d_max＝0.75(σ_L/ρ_L)^0.6ε^-0.4 (75)

式中，k_G为气测传质系数，mol/(Pa·m²·s)；k_L为气测传质系数，m/s；k_A为基于反应速率r_A的一级本征反应速率常数，mol·m³/s；D_G为气相扩散系数；D_L为液相扩散系数；d₃₂为反应器内气泡Sauter平均直径，m；d_min、d_max分别为反应器中最大气泡和最小气泡直径，m；v₃₂为与d₃₂对应的反应器内气泡平均上升速度，m/s；v_G为反应器内的表观气速,m/s；v_L为反应器内的表观液速,m/s；v₀：无限大静止液相中单个气泡上升速度，m/s；S₀：反应器横截面积，m²；D₀：反应器直径，m；Q_L为反应器内液体循环体积流量；L_mix为气泡破碎区长度；D₁为气泡破碎管直径；λ₁为气液体积流量之比，λ₁＝Q_G/Q_L；Q_G为反应器内通气体积流量；P₀为气泡破碎器入口处液体的静压；P₁为气泡破碎区出口气液混合物压力；ε为能量耗散率；μ_L为液体动力粘度；σ_L为液体表面张力；ρ_L为液体密度。实施例3选用的模型主要考虑L_mix小于L_b的情况，因为相反的情况不常见，比较极端。反应器结构参数还需满足：λ₁＝0.1～0.5、K₁＝0.5、L_b＝13D₁。

对于二氧化碳和水体系，当操作条件为：Q_L＝2000L/h(5.56×10^-4m³/s)，气体流量Q_G＝0.2Q_L，T＝298K，P_G0＝1atm；而此体系中液相的物性参数为：ρ_L＝1000kg/m³，μ_L＝8.9×10^-4Pa·s，σ_L＝7.197×10^-4N/m；反应器气泡破碎管直径D₁＝0.02m；E₀表示体系输入的能量，即循环泵铭牌上的额定功率，取E₀＝1000W。根据操作条件和上述模型应用MIR计算得到气泡Sauter平均直径d₃₂＝0.1mm时的反应速率为r_A＝1.42×10^-1mol·m³/s。

表1和表2是在MIR中相同体系的各个参数的计算情况：

表1在MIR反应器中模型公式计算的参数

表2在MIR反应器中模型公式计算的三种阻力(气膜、液膜、本征)

实施例3：

本实施例以图1所示的反应器为例，说明实施例1所述建模方法构建的模型在空气-水体系反应器中的应用，与现有装置反应体系中反应速率r_A相比的优越之处。

表3和表4是相同体系不同粒径下的各个参数的对比情况：

表3不同粒径下模型公式计算的参数

表4不同粒径下模型公式计算的三种阻力(气膜、液膜、本征)

如表3、表4所示，在相同条件下，对于传统反应器产生的直径最小为1mm的气泡相比，MIR产生的气泡直径仅为0.1mm，为原来的1/10，其宏观速率要比前者大了约4倍，而反应阻力逐渐由液膜控制过渡到由本征反应阻力控制。可见气泡直径的减小确实减小了传质阻力和强化了气液传质速率。

Claims (6)

1.一种微界面强化反应器反应速率构效调控模型建模方法，其特征在于，包括：

获取微界面强化反应器反应体系工艺条件，包括催化剂、温度、压力、物料配比；计算气液反应反应速率；

基于Levenspiel的理论，多相体系的反应速率由下式表示：

其中，r_A为反应物A的反应速率；H_A为反应物A的亨利常数；k_G、k_L、k_S分别为气侧、液侧及液固传质系数；a、a_S分别为气液相界面积、液固相界面积；k_A为基于反应速率r_A的一级本征反应速率常数；

为催化剂颗粒内部溶剂B的平均浓度；χ_A为反应速率常数k_A时的一级本征反应有效因子，表征因催化剂孔扩散导致的反应速率降低的程度；f_s为催化剂装载率；P_G为气泡内压；

选定χ_A＝1；

是催化剂表面溶剂的浓度，在数量上与f_s形成反向对应关系，即有：

忽略式(1)分母第三项，将气液反应反应速率方程简化为：

P_G表示为：

式中，ρ_L为液体密度；H₀为反应器内初始液位高度；d₃₂为反应器内气泡Sauter平均直径；v₃₂为与d₃₂对应的反应器内气泡平均上升速度；v_G为反应器内的表观气速；σ_L为液体表面张力；

将(4)带入(3)，得到式(5)：

2.根据权利要求1所述的方法，其特征在于，所述传质系数采用如下方式获取：

气侧传质系数k_G：根据膜理论，气侧传质系数k_G大小由气相扩散系数D_G及有效气膜厚度δ_G决定；

基于Chapman-Enskog动理论，D_G的一般形式为：

式中，T为气体温度；M_A、M_B分别为气体A及溶剂B的摩尔质量；

采用Gedde静止球模型建立k_G理论计算模型：

式中，d₃₂为反应器内气泡Sauter平均直径；t₃₂为直径为d₃₂的气泡在体系中的停留时间，其表达式为：

式中，H₀为反应器内初始液位高度；v₃₂为与d₃₂对应的反应器内气泡平均上升速度；v_G为反应器内的表观气速；

将式(8)代入上式(7)后得：

液侧传质系数k_L：建立基于表面更新理论的渗透模型，具体如下：

式中t_c为液体微元暴露时间；D_L为液相扩散系数；

依据速度滑移模型的定义，将t_c利用式(11)定义：

t_c＝d₃₂/v_S (11)

式中v_s为气泡滑移速度；

将式(11)代入式(10)得：

根据气泡滑移速度v_s的定义：

v_S＝v₃₂-u (13)

式中u为气泡间流体速度；

当液体循环流量为Q_L，对于均匀上升气流，u的表达式为：

式中，v_L为表观液速，

S₀为反应器横截面积，D₀为反应器直径；φ_G为反应器内气含率，假定为φ_G＝v_G/v₃₂；

将式(14)代入式(13)得：

将式(15)代入式(12)得：

式中，液相扩散系数D_L基于修正的Wilke–Chang方程获取，如式(17)所示：

3.根据权利要求2所述的方法，其特征在于，所述气液相界面积a采用如下方式获取：

步骤1：微界面强化反应器气液相界面积通用表达式建立：

假设微界面强化反应器内气泡为刚性球体，设定气泡最大和最小直径分别为d_max和d_min；当反应器稳态工作时，有：

式中，a为反应器内气液相界面积，即单位体积气液混合物中气泡的表面积；x、N、F(x)分别为气泡直径、稳态时反应器内气泡总数及气泡粒径分布概率密度函数；V_T为气液混合物总体积；

设反应器内初始液位高度和反应器直径分别为H₀和D₀，体系气含率为φ_G，计算V_T和反应器横截面积S₀：

将式(19)代入式(18)，得：

依据φ_G的定义：

将式(19)代入式(22)，得：

在稳态情况下，将慢反应体系中气泡单程上升时由于反应消耗的体积及数量忽略不计，则进出反应体系的气体分布基本相等，且气泡逃逸概率f(x)为：

f(x)NF(x)＝n_II(x) (24)

上式中，n_I和I(x)分别为单位时间进入反应器的气泡数量和粒径分布概率密度函数；当反应体系的通气量为Q_G，依据气体体积守恒原则，有：

将式(25)代入(24)并积分，得：

将式(26)代入式(23)后，得：

式(27)中，v_G为表观气速，由下式表示：

v_G＝Q_G/S₀ (28)

假设气泡逸出概率f(x)等于气泡在体系中停留时间的倒数，则有：

f(x)＝v(x)(1-φ_G)/H₀ (29)

式中v(x)是直径为x的气泡上升速度；将式(29)分别代入式(26)、(27)，得：

将式(30)代入式(21)后，得：

假设气泡进入反应器后一般不发生聚并和破裂现象；且由于是慢反应过程，气泡从反应器底部上升至顶部的单程运动过程中，所反应消耗的量可忽略不计，即稳态时反应器内各处的气泡粒径分布与初始气泡分布一致，也就是：

F(x)＝I(x) (32)

将式(32)代入式(31)并进行离散化处理后得到a的计算公式为：

步骤2：微界面强化反应器气含率通用表达式建立：

基于已有a的另一通用表达式，即式(34)：

a＝6φ_G/d₃₂ (34)

式中，d₃₂为气泡Sauter平均直径；对比(33)和(34)，当φ_G定义为下式(35)时，上述两式等价：

式(35)即为φ_G通用计算公式；

步骤3：获取微界面强化反应器气液体系内气泡上升速度v₃₂；

φ_G是体系宏观特征参数，与d₃₂类似，假定：

φ_G＝v_G/v₃₂ (36)

可得：

对于微米级气泡，考虑升力、气泡周围压力梯度、气泡对其周围湍流场响应时间τ_b(s)对曳力的影响，微气泡运动的矢量方程可表示为式(38)：

式中，u为气泡间流体速度；y(t)表示t时刻气泡的位置；g为重力加速度；ω为气泡涡量；τ_b为气泡响应时间；

对于微界面强化反应器的微气泡体系，反应器内为无旋势流场，因此ω＝0；假设气泡进入体系后宏观上为匀速上升，并且稳态时，反应器内湍流各向同性，气泡稳定上升，则Du/Dt＝0；dv₃₂/dt＝0；气泡响应时间τ_b表示为：

式中，v₀为无限大静止液相中单个气泡上升速度；

仅考虑反应器内气泡一维运动情况并结合以上假设，将式(38)简化：

v₃₂＝v₀±u (40)

分析反应器内气泡周围液体的三种运动情况：

a.气泡周围液体宏观竖直向下运动；

此时：

v₃₂＝v₀-u (41)

当体系中气泡较小时，基于反应器工业生产需求，须满足：

v₀≥u (42)

结合式(14)(41)，有：

且：

对式(43)化简后得：

v₃₂ ²+(v_L-v_G-v₀)v₃₂+v_Gv₀＝0 (45)

求解式(45)，获得唯一实解：

b.气泡周围液体宏观运动速度为零；

根据微界面强化反应器中微气泡主导的反应体系，将v₀表达为：

式中，σ_L为液体表面张力；μ_L为液体动力粘度；ρ_L为液体密度；

c.气泡周围液体宏观竖直向上运动：

此时：

v₃₂＝v₀+u (48)

并且：

v₃₂ ²-(v_L+v_G+v₀)v₃₂+v_Gv₀＝0 (49)

整理后得：

根据所述a、b、c，气泡上升速度v₃₂依据液体宏观流动方向表示为：

式中，v_L前的“+”表示气泡周围液体宏观竖直向上运动的情况，“-”则表示气泡周围液体宏观竖直向下运动的情况；上式中v₀依据式(47)计算。

4.根据权利要求1-3任一项所述的方法，其特征在于，所述反应器内气泡Sauter平均直径d₃₂采用如下方式获取：

(1)以微界面强化反应器最大气泡直径d_max和最小气泡直径d_min为自变量，气泡Sauter平均直径d₃₂为因变量，建立d_max、d_min和d₃₂间的关系；具体步骤如下：

设x,m,n分别为反应器气液体系中的气泡粒径、气泡粒径几何自然对数的均值和标准差，获取气泡粒径x的概率密度函数：

气泡粒径满足此分布时的气泡Sauter平均直径d₃₂为：

d₃₂＝exp(m+2.5n²) (53)

气泡粒径x呈对数正态分布，因此lnx的数学期望为：

根据气泡粒径x的概率密度函数绘制气泡粒径概率密度图，当

时，概率密度最大；即此处的一阶导数为0：

对方程(52)求导得到方程(55)：

由(54)、(55)可得：

由于：

将方程(52)代入(57)并化简后可得：

令：

则上式简化为：

方程(59)左端为误差函数，与标准误差函数的差别在于积分限的不同，将式(56)分别代入上述积分上下限，并将方程(59)转化为标准误差函数后可得：

方程(60)中，erf(·)为误差函数；

对于如下形式的误差函数：

采用泰勒级数展开进行近似计算，泰勒级数展开依据误差函数自变量的取值范围不同而采用不同的形式，当z≤4时，erf(z)可展开为：

由于：

当d_max/d_min为1000时：

而根据方程(62)：

因此，当：

即：

时，方程(60)近似成立；

此外，方程(60)成立的条件与n及d_max/d_min的大小有关，且n受d_max/d_min的大小的制约；构建气泡粒径累积概率密度g(n)以考察n和d_max/d_min对方程(60)成立条件的影响，令气泡粒径累积概率密度g(n)为：

绘制g(n)～n关系曲线；获取确保方程(60)成立的n的可取值范围与d_max/d_min的关系；

取不等式(67)的等号条件，即：

由式(56)及(69)确定m和n，进而由方程(53)建立d₃₂基本数学模型；其结果如下：

(2)基于Kolmogorov-Hinze理论，构建微界面强化反应器最大气泡直径d_max、最小气泡直径d_min与反应器参数间的关系；

能使气泡破裂的最小湍流涡尺度是Kolmogorov尺度的11.4～31.4倍，假设此倍率为11.4，由于湍流涡仅能破碎直径大于其尺度的气泡，因此，气泡直径最小值d_min与该湍流涡尺度一致，即：

d_min＝11.4(μ_L/ρ_L)^0.75ε^-0.25 (72)

基于Kolmogorov-Hinze理论，最大气泡直径d_max由下式(73)确定：

d_max＝ε^-2/5(σ_LWe_crit/2ρ_L)^3/5 (73)

其中，ε为能量耗散率；μ_L为液体动力粘度；σ_L为液体表面张力；ρ_L为液体密度；We_crit为气泡破碎临界weber数；

基于气泡破碎的共振理论确定We_crit：

其中，ρ_G为气体密度，g/m³；α₂为气泡体积模量，α₂＝2,3,...；当α₂越大，气泡高阶振动越激烈，气泡就越小，对于超细气泡颗粒选择α₂＝2，即We_crit＝1.24；

此时：

d_max＝0.75(σ_L/ρ_L)^0.6ε^-0.4 (75)

将式(72)、式(75)代入式(71)求解d₃₂。

5.根据权利要求4所述的方法，其特征在于，所述能量耗散率ε采用如下方式获取：

步骤100：将微界面强化反应器总的能量耗散率ε的计算划分为微界面强化反应器内三个不同区域能量耗散率的总和，包括反应器主体区鼓泡区的能量耗散率ε_R，气液破碎区的能量耗散率ε_mix以及气液出口区的能量耗散率ε_pl；

步骤110：其中，反应器主体区鼓泡区的能量耗散率ε_R采用如下方式计算：

反应器气体鼓泡过程中，依据气泡对体系做功，ε_R表示为：

其中，Q_G为反应器内通气体积流量，m³/s；S₀为反应器横截面积，m²；

步骤120：计算气液破碎区的能量耗散率ε_mix：

基于ε_mix经典计算模型，假设气液混合为绝热过程并忽略液体势能变化，忽略气体质量流量，并使能量耗散率的单位统一为W/Kg，获取计算公式如下：

其中，L_mix为气泡破碎区长度，m；P₀、P₁分别为气泡破碎区入口液体静压及出口气液混合物压力，Pa；λ₁为气液体积流量之比：K₁为喷嘴直径与气泡破碎区直径的比值，K₁＝D_N/D₁；S₁为气泡破碎器横截面积，m²；ρ_L为液体密度，kg/m³；Q_L为反应器内液体循环体积流量，m³/s；

λ₁＝Q_G/Q_L (78)

步骤121：计算气泡破碎区入口液体静压P₀及出口气液混合物压力P₁：

忽略气泡破碎区管壁摩擦损失，则：

其中，φ_mix为气泡破碎区气含率，按下式计算：

忽略管道摩擦及喷嘴处能量损失，依据能量守恒原理，体系实际获得的能量E₀为：

即：

由式(79)、(82)得：

步骤122：计算气泡破碎区长度L_mix：

L_mix通过测量气液破碎区管内壁压力突变来确定，或通过如下方式确定：

其中：P_H为气液破碎区上方气压，Pa；ρ_MZ为气液破碎区内气液混合物密度，Kg/m³；v_N为射流口的射流速度，m/s；U_e,max为气液破碎区涡旋最大返回速度，m/s；

P_H由伯努利方程推得：

P_H≈P_G0 (85)

式中，P_G0为供气压力，Pa；

ρ_MZ通过下式计算：

ρ_MZ＝ρ_Gφ_mix+ρ_L(1-φ_mix)≈ρ_L(1-φ_mix) (86)

式中，ρ_G为气体密度，g/m³；

考虑气液破碎区气液混合物流速的影响，U_e,max为射流口射流流速与气液破碎区气液混合物流速的矢量合成的结果，采用下式计算U_e,max：

将式(80)、(85)、(86)、(87)代入式(84)，并化简后可得：

获取反应器气泡破碎器长度L_b，并根据式(88)计算L_mix；

①当L_mix＜L_b时，以式(88)的计算结果为L_mix的实际数值；

②当L_mix≥L_b时，说明射流能量近似完全在气泡破碎器区域耗，则：

L_mix＝L_b (89)

步骤130：计算气液出口区的能量耗散率ε_pl；

假设气泡在气液出口区呈均匀分布状态，此区域的能量耗散速率ε_pl由下式计算：

式中，v_L为表观液速；反应器结构设计时保证λ₁可调范围足够大，通过实验确定反应器基本结构参数间的关系为K₁＝0.5，L_b＝13D₁；代入前述相应表达式并化简可得：

步骤200：确定ε_R、ε_mix和ε_pl各自的数值大小；

步骤210：依据进入反应器的气液流量等于气泡破碎区出口的气液流量平衡原理，得到：

式中，C_L为基于反应器内有效体积πD₀ ²H₀/4的液体循环倍数，即每小时液体循环总体积与反应器有效体积的比值；u₁为气泡破碎器出口气液混合物线速度，m/s；λ₁取值0.1～0.5；

由式(94)可知：

则u₁增大时，反应器横截面积S₀也增大；结合式(76)可知，此时ε_R减小；为对反应器不同区域的能量耗散率进行比较，假设：u₁＝3.0m/s；C_L＝20；H₀＝1.5m；由式(94)可得，当λ₁＝0.1～0.5时：

D₀≈19D₁ (95)

选定D₁数值，计算并比较反应器不同区域在不同喷嘴液速下的能量耗散率，确定与气液破碎区的能量耗散率ε_mix相比，反应器主体区鼓泡区、气液出口区的能量耗散率可忽略不计，即ε_mix≈ε；则整个反应器的能量耗散率ε与反应器结构参数之间的数学关系，可由式(77)计算确定，即：

6.权利要求1-5任一项所述方法在反应器设计中的应用。

Images