本发明属于化工制造、反应器、建模技术领域，具体涉及一种微界面强化反应器传质速率构效调控模型建模方法。
背景技术：
氧化、加氢、氯化等多相反应在化工生产过程中广泛存在，其宏观反应速率一般受制于传质过程。气液反应的传质速率主要受液侧(或气侧)传质系数及气液相界面积a共同影响。已有研究表明，a对体积传质系数的影响程度更大，且容易调控。因此，增大a被视为提高受传质控制的气液反应体系反应效率的特别有效的途径。
气泡sauter平均直径d32是决定a大小的关键参数之一，它们主要受气泡间及气液两相间相互作用力影响。气泡聚并和分裂则分别是上述两种作用力的结果，并影响气泡直径的大小。因此，气泡聚并和破裂作为气泡的介观尺度行为，是决定a大小的深层次原因。关于气泡聚并和分裂行为的研究由来已久，普遍认为能量耗散率和d32是重要的影响因素。事实上，d32能够影响a及体积传质系数大小，是决定气液宏观反应速率的核心因素[1]。研究显示，当d32逐渐减小时，体积传质速率逐渐增大；特别是当d32小于1mm时，体积传质速率随d32的减小以类似于指数形式较快增大。因此，尽可能地减小d32能够强化气液传质并最终增大宏观反应速率。
鼓泡反应器和搅拌-鼓泡反应器是工业上最传统和常用的气液反应器。如px氧化制ta的塔式鼓泡反应器，气泡直径通常大于10mm，乃至几厘米级，其传质界面面积十分有限，因此必须将反应器做得很大，以提高宏观反应速率，同时必须通过增加鼓气量来促进液体湍流，使气含率提高，进而增大界面面积，但此举必然降低空气中氧的利用率，增大压缩机功率和尾气排放，导致能耗过渡和物料损失及环境污染。从湍流动力学角度看，传统上用得最广泛的搅拌-鼓泡式气液反应器内大多形成对气泡宏观运动有影响但对气泡破碎作用甚微的大涡，气泡不能有效破碎，故气泡直径偏大，传质面积受限，以致反应效率偏低。为强化气液传质，塔式鼓泡反应器一般增设气体分布板、静态混合器等内件以加强混合，而搅拌釜则需安装不同结构的搅拌桨或内筒等结构，以增加液层的含气量。尽管如此，这两种反应器内的气泡直径通常为5～20mm，所提供的单位体积中的相界面积均十分有限，一般小于100m2/m3,故反应效率不可能获得突破性提高。因此，工业上经常通过高温高压和加大气量来提高气含率和相界面积，但这对反应过程的能耗、物耗及反应选择性都有重大的负面影响。
由于研发气泡的微破碎技术十分重要，故最近10年来，英、美、德、日等国的大学与研究机构开始关注和研发超细气泡技术[2-11]，但其研究成果具有下列共同缺陷：
(1)采用机械破碎、流体撞击、超声等手段虽可得到一定量微米级尺度的气泡，但气液比(气体体积与液体体积之比)太低，一般低于1％，上限不超过5％。此外，产生微气泡的设备能耗和制造成本也太高。
(2)国内外尚没有基于液相为连续相且高度湍流的微气泡体系特性，提出过系统化的微界面传质强化理论、微气泡测试与表征方法、微界面强化反应器构效调控理论及相关数学模型。
基于上述原因，虽有零星的应用实验结果发表，但尚未有规模化的工业应用报道，尤其在化工制造领域的应用基本上还处于空白状态。
当今时代的化工生产，基于创新、绿色、环保总体考虑，其生存和发展依赖于对材料与过程技术大幅创新。提高反应与分离过程的原子经济性对降低能耗、物耗、增强竞争力至关重要。基于此，我们提出“微界面传质强化反应-精细分离集成系统”新技术，试图从最基本的超细气液颗粒特性研究出发，解决在高度湍流状态下超细颗粒体系中，微界面化学反应器内流体流动、传质、反应、能量转换等构效调控全过程所涉及的理论、技术与应用问题。
本发明涉及的超细气液颗粒指的是超细气泡(或超细液滴)，是颗粒当量直径处在1μm≤d0＜1mm的微米级气液颗粒。在反应体系中，超细气液颗粒形成了超细界面(或微界面)，超细界面的形成大幅强化了传质与反应速率，尤其是受传质控制的反应体系。
经典的气液混合理论一般基于毫米－厘米级气液颗粒特征，目前最为合理的方法是多尺度能量最小化原理(emms)[12]。当前大多研究工作均是针对传统气液反应器内的毫米级上气泡[13,14]，很少涉及超细颗粒体系。针对超细颗粒体系的混合、传质与反应特性，必须建立新的计算模型、测试与表征方法，以及构效调控模型，为此必须研究新的设备结构、能量输入方式与转换模式，从而形成全新的适合于超细颗粒反应体系的计算软件与硬件平台，为我国的过程工业生产技术上一个新台阶提供技术与装备支持。
传统反应器的传质建模是基于实验数据拟合得到体积传质参数(如液相体积传质系数kla、气相体积传质系数kga等)经验或半经验的关联式。此类模型较多，其缺陷在于：所得关联式是在特定条件下(包括特定反应器结构、搅拌方式、气液物性、操作工况、测量方法等)得到的，难以普遍适用。
基于连续介质假设的navier-stokes方程(n-s方程)，理论上可以得到设计参数。但由于工业气液反应体系中气液界面存在复杂的时空变化，即使采用上述方法能获得数值解，其结果的适用性仍令人生疑。主要原因在于求解方程所依赖的边界条件受反应器结构(外壳和内件、搅拌器结构等)和流体运动方式的影响较大，且很难确定。
气液界面传质阻力由kla和kga大小决定。但由于多数气液反应，如加氢和氧化，气体在液相中的溶解度不大，气侧传质阻力较液侧一般要小几个数量级[15,16]，故多数情况下可以忽略气侧传质阻力。因而气液界面传质建模一般是针对反应器结构、气液物性、操作工况和kla之间的关系展开[17-23]。但就传质过程而言，kla是一个宏观参数，它并没有充分反映气液传质界面及其两侧的详细信息。因此，只有将kla中的两项(kl和a)分别加以研究才能获取更多微观信息。为建立普遍适用的气液界面传质模型，还须考虑气侧传质阻力，一概忽略该阻力是不科学的，应依据不同体系的计算或测定结果决定是否可忽略。因此，为建立精确普适性的反应体系气液传质模型，需要对kl、kg和气液相界面积a分别进行研究[24-26]。
经典kg数学模型一般膜理论。据此，决定其大小的参数有温度t、亨利系数ha、气相扩散系数dg及有效气膜厚度δg。由于影响δg大小的因素较多，且不易测量，故尚未见其理论计算模型。除此之外，实际可采用gedde静止球模型[27]：
相对于kg而言，kl的理论研究较多。自上世纪20年代发展至今，已有包括双膜理论[28]、渗透理论[29]及表面更新理论[30]三种经典理论及在它们基础上的一些改进型模型[31-35]，如速度滑移模型[36]，界面湍流理论[37,38]，单涡模型[39]，表面更新伸展模型[40]等。而基于膜理论的kl模型中，有效液膜厚度δl的准确测量存在与δg相同的测量难度的问题，因此用于实际理论计算的意义不大[41,42]。而基于表面更新理论的kl渗透理论模型，则由于液相扩散系数dg和液体微元暴露时间tc(s)可以分别建立相应理论计算模型，故被广泛采用。
所谓构效调控数学模型，是指将超细气液颗粒反应体系之反应效率(能效和物效)与体系理化特性、微界面特性、传质特性和反应器结构用数学方法关联起来，从而实现可通过调整结构参数和操作参数以获得反应过程能效物效的最大化目标，或者在给定反应目标(任务)和能耗物耗下，设计出高效的反应器结构。而对于微界面强化反应器而言，此方面工作在国际上仍为空白。
参考文献
[1]levenspielo.chemicalreactionengineering[m].wileynewyorketc.,1972.
[2]xujh,lisw,chengg,luog..formationofmonodispersemicrobubblesinamicrofluidicdevice[j].aichejournal,2006,52(6):2254-2259.
[3]lipandtsugeh.ozonetransferinanewgas-inducedcontactorwithmicrobubbles[j].journalofchemicalengineeringofjapan,2006,39(11):1213-1220.
[4]muroyamak,imaik,okay,hayashij,masstransferpropertiesinabubblecolumnassociatedwithmicro-bubbledispersions[j].chemicalengineeringscience,201,100:464-473.
[5]maeday,hosokawas,babay,tomiyamaakio.generationmechanismofmicro-bubblesinapressurizeddissolutionmethod[j].experimentalthermalandfluidscience,2015,60:201-207.
[6]hasegawah,nagasakay,kataokah.electricalpotentialofmicrobubblegeneratedbyshearflowinpipewithslits.fluiddynamicsresearch,2008,40(7-8):554-564.
[7]weberjandagblevorf.microbubblefermentationoftrichodermareeseiforcellulaseproduction[j].processbiochemistry,2005,40(2):669-676.
[8]rehmanf,medleygj,bandulasenah,zimmermanwb.fluidicoscillator-mediatedmicrobubblegenerationtoprovidecosteffectivemasstransferandmixingefficiencytothewastewatertreatmentplants[j].environmentalresearch,2015,137:32-39.
[9]strideeandedirisinghem.novelmicrobubblepreparationtechnologies[j].softmatter,2008,4(12):2350.
[10]druzinecd,salzigde,kraumem,czermakp.micro-bubbleaerationinturbulentstirredbioreactors:coalescencebehaviorinpluronicf68containingcellculturemedia[j].chemicalengineeringscience,2015,126:160-168.
[11]李宝璋,尚龙安,姜信真.倒喷射式环流反应器的研究[j].西北大学学报(自然科学版).1989,04:65-69.
[12]chenjh,yangn,gew,lijh.stability-drivenstructureevolution:exploringtheintrinsicsimilaritybetweengas-solidandgas-liquidsystems[j].chinesejournalofchemicalengineering.2012,20(1):167-177.
[13]hinzejo.fundamentalsofthehydrodynamicmechanismofsplittingindispersionprocesses[j].aichejournal.1955,1(3):289-295.
[14]zhongs,zoux,zhangzb,tianhz.aflexibleimageanalysismethodformeasuringbubbleparameters[j].chemicalengineeringscience,2016,141(17):143-153.
[15]akitak,yoshidaf.bubblesize,interfacialarea,andliquid-phasemass-transfercoefficientinbubblecolumns[j].industrial&engineeringchemistryprocessdesignanddvelopment.1974,13(1):84-91.
[16]schumpea,deckwerwd.organicliquidsinabubblecolumn:holdupsandmasstransfercoefficients[j].aichejournal.1987,33(9):1473-1480.
[17]jordanu,schumpea.thegasdensityeffectonmasstransferinbubblecolumnswithorganicliquids[j].chemicalengineeringscience.2001,56(21):6267-6272.
[18]shahyt,kelkarbg,godbolesp,etal.designparametersestimationsforbubblecolumnreactors[j].aichejournal.1982,28(3):353-379.
[19]sadae,kumazawah,leec,etal.gasholdupandmass-transfercharacteristicsinathree-phasebubblecolumn[j].industrial&engineeringchemistryprocessdesignanddevelopment.1986,25(2):472-476.
[20]mandala,kundug,mukherjeed.interfacialareaandliquid‐sidevolumetricmasstransfercoefficientinadownflowbubblecolumn[j].thecanadianjournalofchemicalengineering.2003,81(2):212-219.
[21]kawasey,halardb,moo-youngm.theoreticalpredictionofvolumetricmasstransfercoefficientsinbubblecolumnsfornewtonianandnon-newtonianfluids[j].chemicalengineeringscience.1987,42(7):1609-1617.
[22]vázquezg,cancelama,riverolc,etal.applicationofthedanckwertsmethodinabubblecolumn:effectsofsurfactantsonmasstransfercoefficientandinterfacialarea[j].chemicalengineeringjournal.2000,78(1):13-19.
[23]velanm,ramanujamtk.gas-liquidmasstransferinadownflowjetloopreactor[j].chemicalengineeringscience.1992,47(9):2871-2876.
[24]bouaifim,hebrardg,bastould,etal.acomparativestudyofgashold-up,bubblesize,interfacialareaandmasstransfercoefficientsinstirredgas-liquidreactorsandbubblecolumns[j].chemicalengineeringandprocessing.2001,40(2si):97-111.
[25]kantarcin,borakf,ulgenko.bubblecolumnreactors[j].processbiochemistry.2005,40(7):2263-2283.
[26]martínm,montesfj,galánma.masstransferfromoscillatingbubblesinbubblecolumnreactors[j].chemicalengineeringjournal.2009,151(1-3):79-88.
[27]charpentier,j.-c.,masstransferratesingas-liquidabsorbers.advancesinchemicalengineering,1981.11.
[28]lewiswk,whitmanwg.principlesofgasabsorption[j].industrialandengineeringchemistry.1924,16:1215-1220.
[29]higbier.therateofabsorptionofapuregasintostillliquidduringshortperiodsofexposure[m].1935.
[30]danckwertspv.significanceofliquid-filmcoefficientsingasabsorption[j].industrialandengineeringchemistry.1951,43(6):1460-1467.
[31]kishinevskymk.twoapproachestothetheoreticalanalysisofabsorptionprocesses[j].j.appl.chem.ussr.1955,28:881-886.
[32]kishinevskymk,serebryanskyvt.themechanismofmasstransferatthegas-liquidinterfacewithvigorousstirring[j].j.appl.chem.ussr.1956,29:29-33.
[33]dobbinswe.mechanismofgasabsorptionbyturbulentliquids[m].newyorkuniversity,departmentofcivilengineering,1964.
[34]perlmutterdd.surface-renewalmodelsinmasstransfer[j].chemicalengineeringscience.1961,16(3):287-296.
[35]kulkarniaa.masstransferinbubblecolumnreactors:effectofbubblesizedistribution[j].industrial&engineeringchemistryresearch.2007,46(7):2205-2211.
[36]calderba.ph.reviewseriesno3-gasabsorptionfrombubbles[j].transactionsoftheinstitutionofchemicalengineersandthechemicalengineer.1967,45(8):c209.
[37]levichvg,technicas.physicochemicalhydrodynamics[m].prentice-hallenglewoodcliffs,nj,1962.
[38]daviesjt.turbulentphenomena[j].academic,newyork.1972.
[39]luks,leeyh.masstransferineddiesclosetoair‐waterinterface[j].aichejournal.1986,32(9):1546-1554.
[40]jajueeb,margaritisa,karamanevd,etal.applicationofsurface-renewal-stretchmodelforinterfacemasstransfer[j].chemicalengineeringscience.2006,61(12):3917-3929.
[41]vasconcelosj,rodriguesj,orvalhos,etal.effectofcontaminantsonmasstransfercoefficientsinbubblecolumnandairliftcontactors[j].chemicalengineeringscience.2003,58(8):1431-1440.
[42]colombetd,legendred,cockxa,etal.experimentalstudyofmasstransferinadensebubbleswarm[j].chemicalengineeringscience.2011,66(14):3432-3440.
技术实现要素：
本发明的目的在于，克服现有技术的缺陷，提供微界面强化反应器传质速率构效调控模型建模方法。
为实现上述技术目的，本发明采用如下技术方案：
微界面强化反应器微界面强化反应器一种微界面强化反应器传质速率构效调控模型建模方法，包括：
步骤10：建立气侧传质系数kg计算模型；
根据膜理论，气侧传质系数kg大小由气相扩散系数dg及有效气膜厚度δg决定；
基于chapman-enskog动理论，dg的一般形式为：
式中，t为气体温度；ma、mb分别为气体a及溶剂b的摩尔质量；pg为反应器内气泡内气体平均压力；
采用gedde静止球模型建立kg理论计算模型：
式中，d32为反应器内气泡sauter平均直径；t32为直径为d32的气泡在体系中的停留时间，其表达式为：
式中，h0为反应器内初始液位高度；v32为与d32对应的反应器内气泡平均上升速度；vg为反应器内的表观气速，
将式(3)代入上式(2)后得：
步骤20：建立液侧传质系数计算模型；
建立基于表面更新理论的渗透模型，具体如下：
式中tc为液体微元暴露时间；dl为液相扩散系数；
依据速度滑移模型的定义，将tc利用式(6)定义：
tc＝d32/vs(6)
式中vs为气泡滑移速度；
将式(6)代入式(5)得：
根据气泡滑移速度vs的定义：
vs＝v32-u(8)
式中u为气泡间流体速度；
当液体循环流量为ql，对于均匀上升气流，u的表达式为：
式中，vl为表观液速，
将式(9)代入式(8)得：
将式(10)代入式(7)得：
式中，液相扩散系数dl基于修正的wilke–chang方程获取，如式(12)所示：
基于式(4)，式(11)，式(12)计算获取微界面强化反应器传质速率。
其中，所述反应器内气泡sauter平均直径d32采用如下方式获取：
(1)以微界面强化反应器最大气泡直径dmax和最小气泡直径dmin为自变量，气泡sauter平均直径d32为因变量，建立dmax、dmin和d32间的关系；具体步骤如下：
设x,m,n分别为反应器气液体系中的气泡粒径、气泡粒径几何自然对数的均值和标准差，获取气泡粒径x的概率密度函数：
气泡粒径满足此分布时的气泡sauter平均直径d32为：
d32＝exp(m+2.5n2)(14)
气泡粒径x呈对数正态分布，因此lnx的数学期望(算术平均值)为：
将方程(15)代入(13)得到方程(16)：
由(15)、(16)可得：
由于：
将方程(13)代入(18)并化简后可得：
令：
方程(20)左端为误差函数，与标准误差函数的差别在于积分限的不同，将式(17)分别代入上述积分上下限，并将方程(20)转化为标准误差函数后可得：
方程(21)中，erf(·)为误差函数；
对于如下形式的误差函数：
采用泰勒级数展开进行近似计算，泰勒级数展开依据误差函数自变量的取值范围不同而采用不同的形式，当z≤4时，erf(z)可展开为：
由于：
当dmax/dmin为1000时：
而根据方程(23)：
因此，当：
即：
时，方程(21)近似成立；
此外，方程(21)成立的条件与n及dmax/dmin的大小有关，且n受dmax/dmin的大小的制约；构建气泡粒径累积概率密度g(n)以考察n和dmax/dmin对方程(21)成立条件的影响，令气泡粒径累积概率密度g(n)为：
绘制g(n)～n关系曲线；获取确保方程(21)成立的n的可取值范围与dmax/dmin的关系；
取不等式(28)的等号条件，即：
由式(17)及(30)确定m和n，进而由方程(14)建立d32基本数学模型；其结果如下：
(2)基于kolmogorov-hinze理论，构建微界面强化反应器最大气泡直径dmax、最小气泡直径dmin与反应器参数间的关系；
能使气泡破裂的最小湍流涡尺度是kolmogorov尺度的11.4～31.4倍，假设此倍率为11.4，由于湍流涡仅能破碎直径大于其尺度的气泡，因此，气泡直径最小值dmin与该湍流涡尺度一致，即：
dmin＝11.4(μl/ρl)0.75ε-0.25(33)
基于kolmogorov-hinze理论，最大气泡直径dmax由下式(34)确定：
dmax＝ε-2/5(σlwecrit/2ρl)3/5(34)
其中，ε为能量耗散率；μl为液体动力粘度；σl为液体表面张力；ρl为液体密度；wecrit为气泡破碎临界weber数；
基于气泡破碎的共振理论确定wecrit：
其中，α2为气泡体积模量，α2＝2,3,...；当α2越大，气泡高阶振动越激烈，气泡就越小，对于超细气泡颗粒选择α2＝2，即wecrit＝1.24；
此时：
dmax＝0.75(σl/ρl)0.6ε-0.4(36)
将式(33)、式(36)代入式(32)求解d32。
进一步的，所述能量耗散率ε采用如下方式获取：
步骤100：将微界面强化反应器总的能量耗散率ε的计算划分为微界面强化反应器内三个不同区域能量耗散率的总和，包括反应器主体区鼓泡区的能量耗散率εr，气液破碎区的εmix以及气液出口区的εpl；
步骤110：其中，反应器主体区鼓泡区的能量耗散率εr采用如下方式计算：
反应器气体鼓泡过程中，依据气泡对体系做功，εr表示为：
其中，qg为反应器内通气体积流量，m3/s；s0为反应器横截面积，m2；
步骤120：计算气液破碎区的εmix：
基于εmix经典计算模型，假设气液混合为绝热过程并忽略液体势能变化，忽略气体质量流量，并使能量耗散率的单位统一为w/kg，获取计算公式如下：
其中，lmix为气泡破碎区长度，m；p0、p1分别为气泡破碎区入口液体静压及出口气液混合物压力，pa；λ1为气液体积流量之比：k1为喷嘴直径与气泡破碎区直径的比值，k1＝dn/d1；s1为气泡破碎器横截面积，m2；ρl为液体密度，kg/m3；ql为反应器内液体循环体积流量，m3/s；
λ1＝qg/ql(39)
步骤121：计算气泡破碎区入口液体静压p0及出口气液混合物压力p1：
忽略气泡破碎区管壁摩擦损失，则：
其中，φmix为气泡破碎区气含率，按下式计算：
忽略管道摩擦及喷嘴处能量损失，依据能量守恒原理，体系实际获得的能量e0为：
即：
由式(40)、(43)得：
步骤122：计算气泡破碎区长度lmix：
lmix通过测量气液破碎区管内壁压力突变来确定，或通过如下方式确定：
其中：ph为气液破碎区上方气压，pa；ρmz为气液破碎区内气液混合物密度，kg/m3；vn为射流口的射流速度，m/s；ue,max为气液破碎区涡旋最大返回速度，m/s；
ph由伯努利方程推得：
ph≈pg0(46)
式中，pg0为供气压力，pa；
ρmz通过下式计算：
ρmz＝ρgφmix+ρl(1-φmix)≈ρl(1-φmix)(47)
式中，ρg为气体密度，g/m3；
考虑气液破碎区气液混合物流速的影响，ue,max为射流口射流流速与气液破碎区气液混合物流速的矢量合成的结果，采用下式计算ue,max：
将式(46)(48)代入式(45)，并化简后可得：
获取反应器气泡破碎器长度lb，并根据式(49)计算lmix；
①当lmix＜lb时，以式(49)的计算结果为lmix的实际数值；
②当lmix≥lb时，说明射流能量近似完全在气泡破碎器区域耗，则：
lmix＝lb(50)
步骤130：计算气液出口区的εpl；
假设气泡在气液出口区呈均匀分布状态，此区域的能量耗散速率εpl由下式计算：
反应器结构设计时保证λ1可调范围足够大，通过实验确定反应器基本结构参数间的关系为k1＝0.5，lb＝13d1；代入前述相应表达式并化简可得：
步骤200：确定εr、εmix和εpl各自的数值大小；
步骤210：依据进入反应器的气液流量等于气泡破碎区出口的气液流量平衡原理，得到：
式中，cl为基于反应器内有效体积πd02h0/4的液体循环倍数，即每小时液体循环总体积与反应器有效体积的比值；u1为气泡破碎器出口气液混合物线速度，m/s；λ1取值0.1～0.5；
由式(55)可知：
d0≈19d1(56)
选定d1数值，计算并比较反应器不同区域在不同喷嘴液速下的能量耗散率，确定与气液破碎区的能量耗散率εmix相比，反应器主体区、气液出口区的能量耗散率可忽略不计，即εmix≈ε；则整个反应器的能量耗散率ε与反应器结构参数之间的数学关系，可由式(38)计算确定，即：
进一步的，所述d32对应的反应器内气泡平均上升速度v32采用如下方式获取：
对于微界面强化反应器中的微米级气泡，考虑升力、气泡周围压力梯度、气泡对其周围湍流场响应时间τb(s)对曳力的影响，微气泡运动的矢量方程可表示为式(58)：
式中，u为气泡间流体速度；y(t)表示t时刻气泡的位置；g为重力加速度；ω为气泡涡量；τb为气泡响应时间；
对于微界面强化反应器的微气泡体系，反应器内为无旋势流场，因此ω＝0；假设气泡进入体系后宏观上为匀速上升，并且稳态时，反应器内湍流各向同性，气泡稳定上升，则du/dt＝0；dv32/dt＝0；气泡响应时间τb表示为：
式中，v0为无限大静止液相中单个气泡上升速度；
仅考虑反应器内气泡一维运动情况并结合以上假设，将式(58)简化：
v32＝v0±u(60)
分析反应器内气泡周围液体的三种运动情况：
a.气泡周围液体宏观竖直向下运动；
此时：
v32＝v0-u(61)
当体系中气泡较小时，基于反应器工业生产需求，须满足：
v0≥u(62)
u采用式(9)表示；
结合式(9)(61)，有：
且：
对式(63)化简后得：
v322+(vl-vg-v0)v32+vgv0＝0(65)
求解式(65)，获得唯一实解：
b.气泡周围液体宏观运动速度为零；
根据微界面强化反应器中微气泡主导的反应体系，将fan通过实验拟合获取的v0的表达式化简得：
式中，σl为液体表面张力；μl为液体动力粘度；ρl为液体密度；
c.气泡周围液体宏观竖直向上运动：
此时：
v32＝v0+u(68)
并且：
v322-(vl+vg+v0)v32+vgv0＝0(69)
整理后得：
根据所述a、b、c，气泡上升速度v32依据液体宏观流动方向表示为：
式中，vl前的“+”表示气泡周围液体宏观竖直向上运动的情况，“-”则表示气泡周围液体宏观竖直向下运动的情况；上式中v0依据式(67)计算。
本发明的另一目的在于提供上述方法构建的反应器传质速率构效调控模型，所述传质速率构效调控模型如下：
dmin＝11.4(μl/ρl)0.75ε-0.25(33)
dmax＝0.75(σl/ρl)0.6ε-0.4(36)
式中，d32为反应器内气泡sauter平均直径；dmin、dmax分别为反应器中最大气泡和最小气泡直径；v32为与d32对应的反应器内气泡平均上升速度；vg为反应器内的表观气速；vl为反应器内的表观液速；v0：无限大静止液相中单个气泡上升速度；s0：反应器横截面积；d0：反应器直径；ql为反应器内液体循环体积流量；lmix为气泡破碎区长度；d1为气泡破碎管直径；λ1为气液体积流量之比，λ1＝qg/ql；qg为反应器内通气体积流量；p0为气泡破碎器入口处液体的静压；p1为气泡破碎区出口气液混合物压力；ε为能量耗散率；μl为液体动力粘度；σl为液体表面张力；ρl为液体密度。
本发明的又一目的在于提供上述方法在反应器设计中的应用。
本发明的方法适用于微界面强化反应器，其核心在于气泡破碎器。气泡破碎器的原理是高速射流所携带的气体相互撞击进行能量传递，使气泡破碎，其结构参数有lb、d1，详细结构见附图1，除此之外该反应器的其它结构参数有d0、h0，具体反应器结构相关内容已经公布于发明人在先申请的专利cn10618766a中，本发明中不再赘述。
本发明公开了微界面强化反应器(mir)构效调控模型的构建方法，所谓构效调控数学模型，是指将超细气液颗粒反应体系之反应效率(能效和物效)与体系理化特性、微界面特性、传质特性和反应器结构用数学方法关联起来，从而实现可通过调整结构参数和操作参数以获得反应过程能效物效的最大化目标，或者在给定反应目标(任务)和能耗物耗下，设计出高效的反应器结构。而对于微界面强化反应器而言，此方面工作在国际上仍为空白。
采用本发明的方法构建的微界面强化反应器传质速率构效调控数学模型有如下有益效果：
(1)可以直观的看出传质速率和气泡大小的关系，不仅填补了国际上此方面的空白，为研究微界面体系奠定了理论基础，从而也可通过调整结构参数和操作参数以实现获得反应过程能效物效的最大化目标，或者在给定反应目标(任务)和能耗物耗下，设计出高效的反应器结构。
(2)本发明的方法构建了dmax、dmin和d32的直接计算关系，而不再采用实验拟合的方式获取d32的具体数值，极大减少了在反应器中应用时产生的误差；现有的方法构建的d32模型多针对的是鼓泡反应器(bubblecolumn,bc)和鼓泡搅拌釜反应器(bubbling-stirredreactor,bsr)内的气液体系，或射流泵(gas-liquidjetbump,gljb)内的空气-水体系。而对于工业微界面反应器(mir)，则不一定适用，其原因在于：①mir内气泡破碎机制有别于上述反应器；②工业气液反应体系中可能涉及到高粘度液体，而现有技术并没有考虑液体粘度对d32的影响；而本发明的方法适用于工业微界面反应器(mir)，其通用性更好。
(3)本发明的方法还构建了能量耗散率ε的表达式，ε是决定气泡破裂机制的关键参数。而ε的增大，意味着气泡更容易破裂，即可以得到尺度更小的气泡，但是针对微米级尺度的能量耗散率的表达式，目前相关研究还是空白，本发明的方法通过严谨的推导过程实现了微界面强化反应器能量耗散率ε的数学模型的构建，为指导新型反应器的设计奠定理论基础。
附图说明
图1是一种反应器结构示意图，用于说明本发明建模方法在反应器装置中的应用。
其中1-反应器，2-泵前阀，3-循环泵，4-泵后阀，5-液体流量计，6-换热器，7-气泡破碎器，8-测温仪，9-下降管，10-气体阀，11-气体流量计，12-气相入口，13-压力表，14-液位计；d0-反应器直径，h0-反应器内初始液位高度，d1-气泡破碎管直径，lb-气泡破碎区长度。
具体实施方式
实施例1
本发明的方法具体包括如下步骤：
步骤10：建立气侧传质系数kg计算模型；
根据膜理论，气侧传质系数kg大小由气相扩散系数dg及有效气膜厚度δg决定；
基于chapman-enskog动理论，dg的一般形式为：
式中，t为气体温度；ma、mb分别为气体a及溶剂b的摩尔质量，kg/mol；pg为反应器内气泡内气体平均压力，pa；
对于中高温气体，上式准确性较高，一般误差在5～10％以内。由式可知，当气体组成一定时，dg随温度升高而增大，随气体分压的增大而减小。低密度气体的dg几乎与气体组成无关，一般为5×10-6～1×10-5m2·s-1，较液相扩散系数dl(10-10～10-9m2·s-1)高几个数量级。
δg通常约定平行于界面的速度分量与相应外流速度相差1％的位置为其外边界。但影响δg的因素较多，且不易测量，故依据膜理论建立kg理论计算模型较困难。实际可采用gedde静止球模型，其数学形式如下：
上式中，d32为反应器内气泡sauter平均直径；t32为直径为d32的气泡在体系中的停留时间，s。其表达式为：
式中，h0为反应器内初始液位高度；v32为与d32对应的反应器内气泡平均上升速度；vg为反应器内的表观气速，
将式(3)代入上式(2)后得：
步骤20：建立液侧传质系数计算模型；
kl体现溶于液相中的气体分子在气液边界层内的传递速率大小，本质上，经典理论均基于边界层假设，不同之处在于气体分子的传递机理有差异。双膜理论认为在静止边界层内，气体分子是通过分子扩散传递，而渗透理论及表面更新理论则多考虑到边界层外流场的影响，认为是液体微元(或湍流涡)更新边界层的结果。prandtl边界层理论将它们进行了统一。该理论认为，气液边界层实际由静止的粘滞层和对流层组成，质量传递过程由分子扩散和对流传质机理控制。danckwerts也指出，kl由分子扩散速率和对流传质速率两部分组成。但通常由于气液体系中气泡易变形且存在尾流，其边界层流动实际十分复杂，故基于prandtl边界层理论建立kl精确数学模型仍然十分困难，必须进行合理的简化。
本发明中建立基于表面更新理论的渗透模型，具体如下：
式中tc为液体微元暴露时间；dl为液相扩散系数；
依据速度滑移模型的定义，考虑到气泡周围液体的绕流，表面更新可近似认为是气泡表面原有液体被新鲜液体替换的过程。进一步而言，气泡(其运动方向上的投影直径为d32)相对于其周围液相滑移(滑移速度vs,m.s-1)是气泡表面更新的原因。因此将tc利用式(6)定义：
tc＝d32/vs(6)
将式(6)代入式(5)得：
根据气泡滑移速度vs的定义：
vs＝v32-u(8)
式中u为气泡间流体速度；
当液体循环流量为ql，对于均匀上升气流，u的表达式为：
式中，vl为表观液速，
将式(9)代入式(8)得：
将式(10)代入式(7)得：
式中，液相扩散系数dl基于修正的wilke–chang方程获取，如式(12)所示：
基于式(4)，式(11)，式(12)计算获取微界面强化反应器传质速率。
现有d32的具体数值采用实验拟合的方式获取，在反应器中应用时可能产生很大的误差，本发明设计的d32算法构建了dmax、dmin和d32的直接计算关系，极大地减小了在反应器中应用时的误差。
具体的，其计算步骤包括：
(1)以微界面强化反应器最大气泡直径dmax和最小气泡直径dmin为自变量，气泡sauter平均直径d32为因变量，建立dmax、dmin和d32间的关系；具体步骤如下：
设x,m,n分别为反应器气液体系中的气泡粒径、气泡粒径几何自然对数的均值和标准差，获取气泡粒径x的概率密度函数：
气泡粒径满足此分布时的气泡sauter平均直径d32为：
d32＝exp(m+2.5n2)(14)
气泡粒径x呈对数正态分布，因此lnx的数学期望(算术平均值)为：
将方程(15)代入(13)得到方程(16)：
由(15)、(16)可得：
由于：
将方程(13)代入(18)并化简后可得：
令：
方程(20)左端为误差函数，与标准误差函数的差别在于积分限的不同，将式(17)分别代入上述积分上下限，并将方程(20)转化为标准误差函数后可得：
方程(21)中，erf(·)为误差函数；
对于如下形式的误差函数：
采用泰勒级数展开进行近似计算，泰勒级数展开依据误差函数自变量的取值范围不同而采用不同的形式，当z≤4时，erf(z)可展开为：
由于：
当dmax/dmin为1000时：
而根据方程(23)：
因此，当：
即：
时，方程(21)近似成立；
此外，方程(21)成立的条件与n及dmax/dmin的大小有关，且n受dmax/dmin的大小的制约；构建气泡粒径累积概率密度g(n)以考察n和dmax/dmin对方程(21)成立条件的影响，令气泡粒径累积概率密度g(n)为：
绘制g(n)～n关系曲线；获取确保方程(21)成立的n的可取值范围与dmax/dmin的关系；
取不等式(28)的等号条件，即：
由式(17)及(30)确定m和n，进而由方程(14)建立d32基本数学模型；其结果如下：
(2)基于kolmogorov-hinze理论，构建微界面强化反应器最大气泡直径dmax、最小气泡直径dmin与反应器参数间的关系；
能使气泡破裂的最小湍流涡尺度是kolmogorov尺度的11.4～31.4倍，假设此倍率为11.4，由于湍流涡仅能破碎直径大于其尺度的气泡，因此，气泡直径最小值dmin与该湍流涡尺度一致，即：
dmin＝11.4(μl/ρl)0.75ε-0.25(33)
基于kolmogorov-hinze理论，最大气泡直径dmax由下式(34)确定：
dmax＝ε-2/5(σlwecrit/2ρl)3/5(34)
其中，ε为能量耗散率；μl为液体动力粘度；σl为液体表面张力；ρl为液体密度；wecrit为气泡破碎临界weber数；
基于气泡破碎的共振理论确定wecrit：
其中，α2为气泡体积模量，α2＝2,3,...；当α2越大，气泡高阶振动越激烈，气泡就越小，对于超细气泡颗粒选择α2＝2，即wecrit＝1.24；
此时：
dmax＝0.75(σl/ρl)0.6ε-0.4(36)
将式(33)、式(36)代入式(32)求解d32。
其中，所述能量耗散率ε是决定气泡破裂机制的关键参数。本发明基于微界面强化反应器，设计了更适用于微气泡体系的ε计算方法。
所述能量耗散率ε采用如下方式获取：
步骤100：将微界面强化反应器总的能量耗散率ε的计算划分为微界面强化反应器内三个不同区域能量耗散率的总和，包括反应器主体区鼓泡区的能量耗散率εr，气液破碎区的εmix以及气液出口区的εpl；
步骤110：其中，反应器主体区鼓泡区的能量耗散率εr采用如下方式计算：
反应器气体鼓泡过程中，依据气泡对体系做功，εr表示为：
其中，qg为反应器内通气体积流量，m3/s；s0为反应器横截面积，m2；
步骤120：计算气液破碎区的εmix：
基于εmix经典计算模型，假设气液混合为绝热过程并忽略液体势能变化，忽略气体质量流量，并使能量耗散率的单位统一为w/kg，获取计算公式如下：
其中，lmix为气泡破碎区长度，m；p0、p1分别为气泡破碎区入口液体静压及出口气液混合物压力，pa；λ1为气液体积流量之比：k1为喷嘴直径与气泡破碎区直径的比值，k1＝dn/d1；s1为气泡破碎器横截面积，m2；ρl为液体密度，kg/m3；ql为反应器内液体循环体积流量，m3/s；
λ1＝qg/ql(39)
步骤121：计算气泡破碎区入口液体静压p0及出口气液混合物压力p1：
忽略气泡破碎区管壁摩擦损失，则：
其中，φmix为气泡破碎区气含率，按下式计算：
忽略管道摩擦及喷嘴处能量损失，依据能量守恒原理，体系实际获得的能量e0为：
即：
由式(40)、(43)得：
步骤122：计算气泡破碎区长度lmix：
lmix通过测量气液破碎区管内壁压力突变来确定，或通过如下方式确定：
其中：ph为气液破碎区上方气压，pa；ρmz为气液破碎区内气液混合物密度，kg/m3；vn为射流口的射流速度，m/s；ue,max为气液破碎区涡旋最大返回速度，m/s；
ph由伯努利方程推得：
ph≈pg0(46)
式中，pg0为供气压力，pa；
ρmz通过下式计算：
ρmz＝ρgφmix+ρl(1-φmix)≈ρl(1-φmix)(47)
式中，ρg为气体密度，g/m3；
考虑气液破碎区气液混合物流速的影响，ue,max为射流口射流流速与气液破碎区气液混合物流速的矢量合成的结果，采用下式计算ue,max：
将式(46)(48)代入式(45)，并化简后可得：
获取反应器气泡破碎器长度lb，并根据式(49)计算lmix；
①当lmix＜lb时，以式(49)的计算结果为lmix的实际数值；
②当lmix≥lb时，说明射流能量近似完全在气泡破碎器区域耗，则：
lmix＝lb(50)
步骤130：计算气液出口区的εpl；
假设气泡在气液出口区呈均匀分布状态，此区域的能量耗散速率εpl由下式计算：
反应器结构设计时保证λ1可调范围足够大，通过实验确定反应器基本结构参数间的关系为k1＝0.5，lb＝13d1；代入前述相应表达式并化简可得：
步骤200：确定εr、εmix和εpl各自的数值大小；
步骤210：依据进入反应器的气液流量等于气泡破碎区出口的气液流量平衡原理，得到：
式中，cl为基于反应器内有效体积πd02h0/4的液体循环倍数，即每小时液体循环总体积与反应器有效体积的比值；u1为气泡破碎器出口气液混合物线速度，m/s；λ1取值0.1～0.5；
由式(55)可知：
d0≈19d1(56)
选定d1数值，计算并比较反应器不同区域在不同喷嘴液速下的能量耗散率，确定与气液破碎区的能量耗散率εmix相比，反应器主体区、气液出口区的能量耗散率可忽略不计，即εmix≈ε；则整个反应器的能量耗散率ε与反应器结构参数之间的数学关系，可由式(38)计算确定，即：
本实施例中，所述d32对应的反应器内气泡平均上升速度v32采用如下方式获取：
对于微界面强化反应器中的微米级气泡，考虑升力、气泡周围压力梯度、气泡对其周围湍流场响应时间τb(s)对曳力的影响，微气泡运动的矢量方程可表示为式(58)：
式中，u为气泡间流体速度；y(t)表示t时刻气泡的位置；g为重力加速度；ω为气泡涡量；τb为气泡响应时间；
对于微界面强化反应器的微气泡体系，反应器内为无旋势流场，因此ω＝0；假设气泡进入体系后宏观上为匀速上升，并且稳态时，反应器内湍流各向同性，气泡稳定上升，则du/dt＝0；dv32/dt＝0；气泡响应时间τb表示为：
式中，v0为无限大静止液相中单个气泡上升速度；
仅考虑反应器内气泡一维运动情况并结合以上假设，将式(58)简化：
v32＝v0±u(60)
分析反应器内气泡周围液体的三种运动情况：
a.气泡周围液体宏观竖直向下运动；
此时：
v32＝v0-u(61)
当体系中气泡较小时，基于反应器工业生产需求，须满足：
v0≥u(62)
u采用式(9)表示；
结合式(9)(61)，有：
且：
对式(63)化简后得：
v322+(vl-vg-v0)v32+vgv0＝0(65)
求解式(65)，获得唯一实解：
b.气泡周围液体宏观运动速度为零；
根据微界面强化反应器中微气泡主导的反应体系，将fan通过实验拟合获取的v0的表达式化简得：
式中，σl为液体表面张力；μl为液体动力粘度；ρl为液体密度；
c.气泡周围液体宏观竖直向上运动：
此时：
v32＝v0+u(68)
并且：
v322-(vl+vg+v0)v32+vgv0＝0(69)
整理后得：
根据所述a、b、c，气泡上升速度v32依据液体宏观流动方向表示为：
式中，vl前的“+”表示气泡周围液体宏观竖直向上运动的情况，“-”则表示气泡周围液体宏观竖直向下运动的情况；上式中v0依据式(67)计算。
实施例2
本实施例以图1所示的反应器为例，说明实施例1所述建模方法构建的模型在二氧化碳和水体系反应器中的应用。图1的反应器结构可为现有反应器的结构，仅采用本发明的方法进行参数设计，本发明中对反应器的结构不再赘述。
根据实施例1构建的传质速率构效调控模型如下：
dmin＝11.4(μl/ρl)0.75ε-0.25(33)
dmax＝0.75(σl/ρl)0.6ε-0.4(36)
式中，d32为反应器内气泡sauter平均直径，m；dmin、dmax分别为反应器中最大气泡和最小气泡直径，m；v32为与d32对应的反应器内气泡平均上升速度，m/s；vg为反应器内的表观气速,m/s；vl为反应器内的表观液速,m/s；v0：无限大静止液相中单个气泡上升速度，m/s；s0：反应器横截面积，m2；d0：反应器直径，m；ql为反应器内液体循环体积流量；lmix为气泡破碎区长度；d1为气泡破碎管直径；λ1为气液体积流量之比，λ1＝qg/ql；qg为反应器内通气体积流量；p0为气泡破碎器入口处液体的静压；p1为气泡破碎区出口气液混合物压力；ε为能量耗散率；μl为液体动力粘度；σl为液体表面张力；ρl为液体密度。
实施例3选用的模型主要考虑lmix小于lb的情况，因为相反的情况不常见，比较极端。反应器结构参数还需满足：λ1＝0.1～0.5、k1＝0.5、lb＝13d1。
对于二氧化碳和水体系，当操作条件为：ql＝2000l/h(5.56×10-4m3/s)，气体流量qg＝0.2ql，t＝298k，pg0＝1atm；而此体系中液相的物性参数为：ρl＝1000kg/m3，μl＝8.9×10-4pa·s，σl＝7.197×10-4n/m；反应器气泡破碎管直径d1＝0.02m；e0表示体系输入的能量，即循环泵铭牌上的额定功率，取e0＝1000w。根据操作条件和上述模型应用mir计算得到气泡sauter平均直径d32＝0.1mm时的气侧和液侧的传质速率分别为kg＝1.78×10-6mol/pa·m3·s，kl＝5.75×10-4m/s。
levenspiel认为，多相体系的宏观反应速率可由下式表示：
经化简后的气液反应宏观速率方程可简化为：
表1和表2是在mir中相同体系的各个参数的计算情况：
表1在mir反应器中模型公式计算的参数
表2在mir反应器中模型公式计算的三种阻力(气膜、液膜、本征)
实施例3：
本实施例以图1所示的反应器为例，说明实施例1所述建模方法构建的模型在空气-水体系反应器中的应用，与现有装置产生的气液传质速率kl和kg相比的结果。
表3和表4是相同体系不同粒径下的各个参数的对比情况：
表3不同粒径下模型公式计算的参数
表4不同粒径下模型公式计算的三种阻力(气膜、液膜、本征)
如表3、表4所示，在相同条件下，通过模型计算，对于传统反应器产生的直径最小为1mm的气泡相比，mir产生的气泡直径仅为0.1mm，为原来的1/10，通过改变了气液传质系数，最终影响了其反应的宏观反应速率，使其增大了将近4倍。此结果在工业生产中具有指导性的意义。