CN107335390B - 微界面强化反应器相界面积构效调控模型建模方法 - Google Patents

Abstract

本发明涉及微界面强化反应器相界面积构效调控模型建模方法，依次构建了微界面强化反应器气液相界面积通用表达式、微界面强化反应器气含率通用表达式和微界面强化反应器气液体系内气泡上升速度表达式，获取微界面强化反应器相界面积的构效调控模型，填补了现有技术的空白。采用本发明的方法可将超细气液颗粒反应体系之反应效率(能效和物效)与体系理化特性、微界面特性、传质特性和反应器结构用数学方法关联起来，从而实现可通过调整结构参数和操作参数以获得反应过程能效物效的最大化目标，或者在给定反应目标(任务)和能耗物耗下，设计出高效的反应器结构。

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]，很少涉及超细颗粒体系。针对超细颗粒体系的混合、传质与反应特性，必须建立新的计算模型、测试与表征方法，以及构效调控模型，为此必须研究新的设备结构、能量输入方式与转换模式，从而形成全新的适合于超细颗粒反应体系的计算软件与硬件平台，为我国的过程工业生产技术上一个新台阶提供技术与装备支持。

传统反应器的传质建模是基于实验数据拟合得到体积传质参数(如液相体积传质系数k_La、气相体积传质系数k_Ga等)经验或半经验的关联式。此类模型较多，其缺陷在于：所得关联式是在特定条件下(包括特定反应器结构、搅拌方式、气液物性、操作工况、测量方法等)得到的，难以普遍适用。

基于连续介质假设的Navier-Stokes方程(N-S方程)，理论上可以得到设计参数。但由于工业气液反应体系中气液界面存在复杂的时空变化，即使采用上述方法能获得数值解，其结果的适用性仍令人生疑。主要原因在于求解方程所依赖的边界条件受反应器结构(外壳和内件、搅拌器结构等)和流体运动方式的影响较大，且很难确定。

还有一种计算相界面积的通用表达式为：

其中φ_G为反应器内气含率；采用上式计算时必须注意两点：其一是上式仅为数理统计意义上的结果，并不是严格数学模型；其二为该式虽然支持气泡粒径分布近似于Maxwell-Boltzmann分布特征的气液鼓泡反应器中a的计算，但对近一二十年发展起来的新型强化反应器以及在特殊搅拌方式下工作的反应体系是否也适用，仍有待研究证实。

气液界面传质阻力由k_La和k_Ga大小决定。但由于多数气液反应，如加氢和氧化，气体在液相中的溶解度不大，气侧传质阻力较液侧一般要小几个数量级^[15,16]，故多数情况下可以忽略气侧传质阻力。因而气液界面传质建模一般是针对反应器结构、气液物性、操作工况和k_La之间的关系展开^[17-23]。但就传质过程而言，k_La是一个宏观参数，它并没有充分反映气液传质界面及其两侧的详细信息。因此，只有将k_La中的两项(k_L和a)分别加以研究才能获取更多微观信息。为建立普遍适用的气液界面传质模型，还须考虑气侧传质阻力，一概忽略该阻力是不科学的，应依据不同体系的计算或测定结果决定是否可忽略。因此，为建立精确普适性的反应体系气液传质模型，需要对k_L、k_G和气液相界面积a分别进行研究^[24-26]。

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

参考文献

[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 F68containingcell 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]Akita K,Yoshida F.Bubble size,interfacial area,and liquid-phasemass-transfer coefficient in bubble columns[J].Industrial&EngineeringChemistry Process Design and Dvelopment.1974,13(1):84-91.

[16]Schumpe A,Deckwer W D.Organic liquids in a bubble column:holdupsand mass transfer coefficients[J].AIChE journal.1987,33(9):1473-1480.

[17]Jordan U,Schumpe A.The gas density effect on mass transfer inbubble columns with organic liquids[J].Chemical Engineering Science.2001,56(21):6267-6272.

[18]Shah Y T,Kelkar B G,Godbole S P,et al.Design parametersestimations for bubble column reactors[J].AIChE Journal.1982,28(3):353-379.

[19]Sada E,Kumazawa H,Lee C,et al.Gas holdup and mass-transfercharacteristics in a three-phase bubble column[J].Industrial&EngineeringChemistry Process Design and Development.1986,25(2):472-476.

[20]Mandal A,Kundu G,Mukherjee D.Interfacial Area and Liquid‐SideVolumetric Mass Transfer Coefficient in a Downflow Bubble Column[J].TheCanadian Journal of Chemical Engineering.2003,81(2):212-219.

[21]Kawase Y,Halard B,Moo-Young M.Theoretical prediction ofvolumetric mass transfer coefficients in bubble columns for Newtonian andnon-Newtonian fluids[J].Chemical Engineering Science.1987,42(7):1609-1617.

[22]Vázquez G,Cancela M A,Riverol C,et al.Application of theDanckwerts method in a bubble column:Effects of surfactants on mass transfercoefficient and interfacial area[J].Chemical Engineering Journal.2000,78(1):13-19.

[23]Velan M,Ramanujam T K.Gas-liquid mass transfer in a down flow jetloop reactor[J].Chemical engineering science.1992,47(9):2871-2876.

[24]Bouaifi M,Hebrard G,Bastoul D,et al.A comparative study of gashold-up,bubble size,interfacial area and mass transfer coefficients instirred gas-liquid reactors and bubble columns[J].Chemical Engineering andProcessing.2001,40(2SI):97-111.

[25]Kantarci N,Borak F,Ulgen K O.Bubble column reactors[J].ProcessBiochemistry.2005,40(7):2263-2283.

[26]Martín M,Montes F J,Galán M A.Mass transfer from oscillatingbubbles in bubble column reactors[J].Chemical Engineering Journal.2009,151(1-3):79-88.

发明内容

本发明的目的在于，克服现有技术的缺陷，提供一种微界面强化反应器相界面积构效调控模型建模方法。

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

一种微界面强化反应器相界面积构效调控模型建模方法，包括：

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

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

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

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

将式(2)代入式(1)，得：

依据φ_G的定义：

将式(2)代入式(5)，得：

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

f(x)NF(x)＝n_II(x) (7)

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

将式(8)代入(7)并积分，得：

将式(9)代入式(6)后，得：

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

v_G＝Q_G/S₀ (11)

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

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

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

将式(13)代入式(4)后，得：

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

F(x)＝I(x) (15)

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

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

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

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

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

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

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

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

φ_G＝v_G/v₃₂ (19)

可得：

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

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

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

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

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

v₃₂＝v₀±u (23)

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

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

此时：

v₃₂＝v₀-u (24)

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

v₀≥u (25)

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

式中，v_L为表观液速，采用下式计算；

结合式(19)(24)(26)，有：

且：

对式(28)化简后得：

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

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

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

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

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

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

此时：

v₃₂＝v₀+u (43)

并且：

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

整理后得：

根据所述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²) (48)

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

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

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

将方程(49)代入(47)得到方程(50)：

由(49)、(50)可得：

由于：

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

令：

则上式简化为：

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

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

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

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

由于：

当d_max/d_min为1000时：

而根据方程(57)：

因此，当：

即：

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

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

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

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

由式(51)及(64)确定m和n，进而由方程(48)建立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 (67)

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

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

其中，ε为能量耗散率；μ_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 (70)

将式(67)、式(70)代入式(66)求解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²) (48)

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

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

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

将方程(49)代入(47)得到方程(50)：

由(49)、(50)可得：

由于：

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

令：

则上式简化为：

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

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

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

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

由于：

当d_max/d_min为1000时：

而根据方程(57)：

因此，当：

即：

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

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

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

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

由式(51)及(64)确定m和n，进而由方程(48)建立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 (67)

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

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

其中，ε为能量耗散率；μ_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 (70)

将式(67)、式(70)代入式(66)求解d₃₂。

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

步骤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 (73)

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

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

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

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

即：

由式(74)、(77)得：

步骤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 (80)

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

ρ_MZ通过下式计算：

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

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

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

将式(80)(82)代入式(79)，并化简后可得：

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

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

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

L_mix＝L_b (84)

步骤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；

由式(89)可知：

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

D₀≈19D₁ (90)

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

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

本发明的另一目的在于提供利用上述建模方法构建的相界面积构效调控模型，所述相界面积构效调控模型如下：

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

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

式中，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为液体密度。

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

由Levenspiel的宏观反应速率方程可知，对于气液反应，气液相界面积a的大小是决定界面传质阻力及宏观反应速率大小的最关键因素，同时亦对反应选择性产生重要影响。因此，a的数学建模与计算是实现气液反应强化和调控的重要途径。

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

(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

本实施例说明本发明的相界面积构效调控模型建模方法。

相界面积a是气液反应体系的关键特征参数之一，有关a的数学建模客观上应针对气泡(液滴)颗粒粒径呈一定离散分布特征的体系。大量实验及理论研究均表明，非均相反应器中气泡(某种液滴)颗粒粒径均呈离散对数正态分布。因此，本文假设这种分布普遍存在于各类非均相反应器中，且存在m类气泡大小。

为便于数学处理，首先假定体系中气泡粒径呈连续分布，由此得到a的连续函数形式，然后进行离散化分析，得到a的通用模型。由于在实际体系中，同时存在于液相中的反应物和产物，使气泡(或某一液相)表面呈污染状态。因此，为简化起见，假设超细气泡(液滴)为刚性球体。

以下基于气液体系进行数学模型推导，液液体系的模型推导过程与之基本相仿。设定气泡最大和最小直径分别为d_max(m)和d_min(m)。依据a的定义，当反应器稳态工作时，有：

上式中，a为反应器内气液相界面积，即单位体积气液混合物中气泡的表面积，m²/m³，x、N、F(x)分别为气泡直径(m)、稳态时反应器内气泡总数及气泡粒径分布概率密度函数。V_T为气液混合物总体积(m³)。若反应器横截面积为S₀(m²)，体系气含率为φ_G，反应器内初始液位高度和反应器直径分别为H₀(m)和D₀(m)，则：

将式(2)代入式(1)，可得

依据φ_G的定义：

将式(2)代入式(5)，可得

在稳态情况下,可以认为慢反应体系中气泡单程上升时由于反应消耗的体积及数量可以忽略不计，故进出反应器体系的气体分布基本相等，且气泡逃逸概率(单位时间内反应器内直径为x的气泡从液面上方逃逸的概率)f(x)为：

f(x)NF(x)＝n_II(x) (7)

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

将式(8)代入(7)并积分，可得

将式(9)代入式(6)后，可得

式(10)中，v_G为表观气速,m/s，由下式表示：

v_G＝Q_G/S₀ (11)

对于均匀气泡流，Lee等认为f(x)为常数的假设并不合理。因为对于实际体系，必定存在一定的气泡粒径分布，而不同大小的气泡其上升速度可能存在较大差异(如毫米级气泡和微米级气泡)。假设气泡逸出概率f(x)等于气泡在体系中停留时间的倒数，则有:

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

上式中v(x)是直径为x的气泡上升速度，m/s。将式(12)分别代入式(9)、(10)，得:

将式(13)代入式(4)后，可得

由于气液反应器内气泡粒径存在上下限，且呈一定离散分布特征，以对数正态分布居多，故须对上式进行离散化处理。同时，虽然v_G、v(x)、I(x)和F(x)决定a的大小，且操作压力^[8]、φ_G对F(x)有影响，但迄今为止后三者与反应器结构参数间的关系仍然缺乏相应的数据。

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

F(x)＝I(x) (15)

将式(15)代入式(14)并进行离散化处理后可得，a的计算公式为：

式(16)即为基于概率论推导得到的气液相界面积a的通用数学表达式。

可以看出，a与v_G、v(x)、F(x)、m四个参数有关。由于气泡粒径存在一定分布特征，故a与v_G一般为非线性关系。

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

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

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

上式中，d₃₂为气泡Sauter平均直径，m。对比(16)和(17)两式可知，当φ_G定义为式(18)时，上述两式等价：

式(18)是φ_G通用计算公式。

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

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

φ_G＝v_G/v₃₂ (19)

其中，v₃₂为与d₃₂对应的体系气泡平均上升速度，m/s。

由以上分析可知：

根据牛顿运动定律，v₃₂由气液两相界面的受力情况决定。后者不仅与气泡大小及其界面性质有关，还受气泡周围液体复杂湍流状况(后者常称之为气泡群效应)的影响。由于湍流具有混沌特性，客观而言，实际反应中气泡运动具有绝对不稳定性。一些研究表明，气泡在体系中停留一段时间后气泡上升速度处于较弱的振荡状态。将v₃₂与气泡曳力系数C_D进行关联，以综合反映各种因素对v₃₂的影响，是研究气泡运动的普遍方法。该法可通过迭代计算较准确地计算v₃₂，但反应器结构参数对v₃₂的影响并不能直观体现。另一种确定v₃₂的思路是，以一定形式的φ_G函数体现各种因素对v₃₂的影响，但多数为经验或半经验形式，难具普遍性是其共同的缺陷。另一方面，当气泡尺度由毫米级到微米级变化时，气泡运动影响因素的重要性将发生改变，具体体现在气泡运动方程具体形式上的差异。对于毫米级气泡，其普遍意义的运动方程可依据Boussinesq-Basset方程，考虑气泡所受浮力、曳力、虚拟质量力、历史力得到^[24]；而对于微米级气泡，除应考虑上述各力外，还应考虑升力、气泡周围压力梯度、气泡对其周围湍流场响应时间τ_b(s)对曳力的影响^[25-29]。

微气泡运动的矢量方程可表示为^[30-34]：

上式中，黑体字母表示矢量。其中，u为气泡间流体速度(隙间流速)，m/s；y(t)表示t时刻气泡的位置；g为重力加速度，m/s²；ω为气泡涡量，s^-1；τ_b为气泡响应时间，s。上式右端第一项为湍流场流体物质导数，体现压力梯度、气泡表面粘性应力、附加质量力的共同影响，后三项分别为曳力、气泡内气体重力、气泡周围湍流涡对气泡运动的影响^[35]。

对于微气泡体系而言，其特征在于：(1)反应器内液体为剧烈湍流，因此，微气泡周围液体的流动可视为势流；(2)微气泡可视为表面被污染的理想硬球。依据上述特征，可对式(21)作如下简化：由于反应器内为无旋势流场，因此,ω＝0；研究表明，气泡进入体系后的上升过程一般不可能为匀速，经历一段时间后其瞬时速度即在一定速度大小附近振荡。此段时间相对于气泡整个停留时间而言较短，因此可认为气泡进入体系后宏观上为匀速上升。同时假设稳态时，反应器内湍流各向同性，气泡稳定上升，因此，Du/Dt＝0；dv₃₂/dt＝0；

气泡响应时间τ_b可表示为：

上式中，v₀为无限大静止液相中单个气泡上升速度，m/s。大量实验表明，对于微气泡，由于其内气体循环对气液界面有一定影响，其运动趋势基本符合颗粒Stokes公式，仅仅是绝对大小上存在一定差异。

当仅考虑反应器内气泡一维运动情况并结合以上假设，式(21)可简化为(以竖直向上方向为正方向)：

v₃₂＝v₀±u (23)

上式中，“+”表示液体宏观竖直向上运动；反之，表示液体宏观竖直向下运动。

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

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

此情况下：

v₃₂＝v₀-u (24)

当体系中气泡较小时，有可能出现v₀<|u|的情况。此时，体系中气泡将整体向反应器下方运动并进入液体外循环管路，由此带来工业生产中的各种后续问题。因此须保证：

v₀≥u (25)

当液体循环流量为Q_L，m³/s,对于均匀上升气泡流，u可表示为：

上式中，v_L为表观液速，m/s：

结合式(19)(24)(26)，有

且：

对式(28)化简后可得：

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

关于v₃₂的方程(30)有解的条件为：

在微气泡主导的气液反应器中，气液表观流速必须满足下列关系：

v_L＞v_G＞0 (33)

因此，式(32)在数学上是不成立的，式(31)即成为方程(30)的唯一有解条件。

根据Stokes定律，微气泡的运动速度应遵从下列方程：

上式中，ρ_L为液体密度，kg/m³；μ_L为液体动力粘度，Pa.s。将式(34)代入式(31)，可得

式(35)表明，体系中平均最小气泡(d₃₂)_min均与μ_L、v_G和v_L的平方根成正比。

此时，方程(30)存在两个实根：

在数学上，式(36)两个实根只有一个是方程(30)的实解。

当反应器内仅有单个气泡时，应有下列关系式成立：

v₃₂＝v₀ (37)

此时方程(36)简化为：

式(38)即为方程(30)的唯一实解，它也是微气泡主导的气液反应体系气泡平均运动速度v₃₂的通用方程。Nicklin表达式(39)仅为本方程的一个特例：

v₃₂＝v_G+v₀-v_L (38)

本研究发现，当v_G足够小时，式(38)即简化为式(39)。后者与前者存在差异的原因在于，Nicklin在推导式(39)时，假设体系为均匀气泡流。实际上，众多的研究已经达成如下共识：对于常见的毫米级气泡体系，v_G是影响气液两相流型的关键因素，当v_G较小时，体系为均匀气泡流，当v_G继续增大时，依次形成非均匀的弹状流、气液团状流及环状流。但对于微米级气泡体系，气泡聚并的概率很小，在较大v_G下仍可能为均匀气泡流。因此，式(39)仅是式(38)的特例。

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

前已述及，气泡在运动过程中内外存在气体循环及诱导湍流等复杂现象，故Fan等通过大量实验数据拟合得到下列v₀的表达式^[43]：

Mo＝gμ_L ⁴/ρ_Lσ_L ³；d_e'＝d_e(ρ_Lg/σ_L)^1/2 (41)

对于微气泡主导的反应体系，(原因)故上式可简化为：

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

此时：

v₃₂＝v₀+u (43)

并且：

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

整理后得：

综上所述，对于超细气泡体系而言，气泡上升速度v₃₂依据液体宏观流动方向可表示为：

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

就气液反应体系而言，无论是搅拌釜式反应器还是塔式鼓泡反应器，其液相出口(液体出料或液相循环)宜设置在反应器上部(液面下方)。因此，在液体宏观向上运动的情况下，单个气泡的运动速度应同时遵从式(29)和式(31)：

气泡Sauter平均直径d₃₂是决定相界面积a大小的关键参数之一，现有d₃₂的具体数值采用实验拟合的方式获取，在反应器中应用时可能产生很大的误差，本发明设计的d₃₂算法构建了d_max、d_min和d₃₂的直接计算关系，极大地减小了在反应器中应用时的误差。

具体的，其计算方法步骤包括：

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

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

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

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

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

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

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

将方程(49)代入(47)得到方程(50)：

由(49)、(50)可得：

由于：

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

令：

则上式简化为：

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

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

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

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

由于：

当d_max/d_min为1000时：

而根据方程(57)：

因此，当：

即：

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

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

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

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

由式(51)及(64)确定m和n，进而由方程(48)建立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 (67)

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

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

其中，ε为能量耗散率；μ_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 (70)

将式(67)、式(70)代入式(66)求解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 (73)

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

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

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

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

即：

由式(74)、(77)得：

步骤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 (80)

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

ρ_MZ通过下式计算：

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

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

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

将式(80)(82)代入式(79)，并化简后可得：

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

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

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

L_mix＝L_b (84)

步骤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；

由式(89)可知：

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

D₀≈19D₁ (90)

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

实施例2

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

构建的相界面积构效调控模型如下：

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

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

式中，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为液体密度。

实施例2选用的模型主要考虑L_mix小于L_b的情况，因为相反的情况不常见，比较极端。

根据实施例1，反应器结构参数还应满足：λ₁＝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时的相界面积为a＝5828.45m²/m³。

Levenspiel认为，多相体系的宏观反应速率可由下式表示：

经化简后的气液反应宏观速率方程可简化为：

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

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

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

实施例3：

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

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

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

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

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

Claims (3)

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

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

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

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

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

将式(2)代入式(1)，得：

依据φ_G的定义：

将式(2)代入式(5)，得：

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

f(x)NF(x)＝n_II(x) (7)

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

将式(8)代入(7)并积分，得：

将式(9)代入式(6)后，得：

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

v_G＝Q_G/S₀ (11)

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

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

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

将式(13)代入式(4)后，得：

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

F(x)＝I(x) (15)

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

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

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

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

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

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

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

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

φ_G＝v_G/v₃₂ (19)

可得：

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

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

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

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

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

v₃₂＝v₀±u (23)

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

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

此时：

v₃₂＝v₀-u (24)

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

v₀≥u (25)

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

式中，v_L为表观液速，采用下式计算；

结合式(19)(24)(26)，有：

且：

对式(28)化简后得：

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

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

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

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

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

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

此时：

v₃₂＝v₀+u (43)

并且：

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

整理后得：

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

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

2.根据权利要求1所述的方法，其特征在于，所述反应器内气泡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²) (48)

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

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

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

将方程(49)代入(47)得到方程(50)：

由(49)、(50)可得：

由于：

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

令：

则上式简化为：

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

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

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

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

由于：

当d_max/d_min为1000时：

而根据方程(57)：

因此，当：

即：

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

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

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

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

由式(51)及(64)确定m和n，进而由方程(48)建立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 (67)

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

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

其中，ε为能量耗散率；μ_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 (70)

将式(67)、式(70)代入式(66)求解d₃₂。

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

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

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

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

其中，v_G为反应器内的表观气速；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 (73)

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

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

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

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

即：

由式(74)、(77)得：

步骤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 (80)

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

ρ_MZ通过下式计算：

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

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

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

将式(80)(82)代入式(79)，并化简后可得：

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

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

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

L_mix＝L_b (84)

步骤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；

由式(89)可知：

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

D₀≈19D₁ (90)

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

Images