微反应器当量直径一般为10~2 000 µm，具有体积小、比表面积大等特点，在生物、化学、电子、机械等领域获得广泛应用^[1-5]. 微反应器管道中的液体流动速度相对较低，不易进行快速混合，大多以无源方法通过改变管道长度或增加湍能强度，以促进液体混合^[6-8]. 在微反应器内部设计分支结构会导致加工成本增加，操作性差；复杂结构或微结构尺寸管道易被固体颗粒阻塞或黏稠物填塞，在管道弯折、交叉等处表现尤为明显. Kashid等^[5-9]的研究表明：固体颗粒的轴径向迁移、团聚现象以及返混行为会在一定程度上导致催化剂失活，引起微反应器管道内部阻塞. 当前关于多相流液-液微反应器、液-固微反应器、气-液-固微反应器等方面的研究，主要涉及流体力学特性、多相流耦合、流型识别等方面，对微反应器特征尺度下的流体流动以及混合特性报道相对较少^[10-11]. 微反应器内流体流动过程中，黏性力起主导作用，液体间的混合需要依赖于分子扩散，混合特性及传质效率较低；因此，在微反应器工作过程中引入外部动力源，改变内部流场分布，是提高工作效能的有效手段.

超声波具有穿透性强、能量密度高、安全可靠等优点，广泛应用于石化、冶金、机械等领域. 鉴于微反应器的几何特性，将超声波激振加载到微反应器中，既可以抑制微反应器阻塞问题，也可以增强混合特性和强化作用. Rahimi等^[12-13]的研究表明，在超声波激振作用下的碘化钾降解率提高约50%，超声波频率、功率等对混合有明显的影响. Parvizian等^[14]研究发现，超声波激振造成的空化特性数值模拟中空化泡破裂产生冲击流和微射流，可以提高局部湍动能与混合强度. Lee等^[15]研究超声波诱导晶体成核的机理发现，频率为20 kHz的超声波能够降低晶体尺寸分布，提高分离率. Yue等^[16]利用超声波激振调控矩形微反应器内部的气液流动特性和传质，仿真结果表明，超声波激振能够有效提高微反应器内部流场的混合强度. Ito等^[17]研究表明，高频超声波能够促进湍流混合和化学反应进程. Balachandran等^[18]研究不同超声波频率对生姜中姜辣素提取率的影响，结果表明当超声波频率为20 kHz时提取率最高. 陈光文等^[19]揭示了多相流动体系中的混合与传质强化机制，与冯浩等^[20]对不同流量速度进行实验观察扰流作用下两相流流动状态相吻合，可以通过改变原混合液层流状态达到湍流效果实现快速混合. Kockmann等^[21]研究T型微反应器在高速流动下的扰流和化学反应，发现入口流速增大管道交汇处对流作用明显雷诺数升高具有一致性；在管道弯折处颗粒团聚现象与钟佳奇^[22]磨粒群分布和动力学特性研究结果相吻合；在超声波促进混合方面和传质等方面得到与Akbari等^[23-24]研究的相似实验现象.

针对上述问题，本文提出超声波耦合作用下的微反应器计算流体力学与离散元耦合（CFD-DEM）建模与调控方法；利用超声波激振技术，强化管道内部混合、传质以及抑制管道内部颗粒阻塞，实现声场与流场、固体相的多物理场耦合，有效进行流场分布调控；结合分形方法，对反应器内部的流场状态进行识别，为流场在线监控与闭环控制提供理论依据与技术支持.

根据多相流混合特性及多场耦合作用效果，以T型汇流微反应器为研究对象，将复杂的流道结构进行化简，如图1所示. 微系统中流体的层流状态是相对于大部分宏观系统中湍流状态的最显著的特点. 基于这一点，微系统中的流体行为更易于采用数值模拟. 相反，湍流状态的流体只能用一些半经验公式进行描述. 在微化工系统中通过合理假设得到的数学模型比宏观系统的应用范围更广. 此外，微化工系统中的流体往往与基于宏观系统得到的连续性假设发生矛盾，这种矛盾在多相流系统中表现得尤为明显，诸如部分学者认为由层流到湍流过渡的临界雷诺数减小、摩擦阻力系数远高于充分发展层流理论计算值；另一些研究结果恰恰相反，认为摩擦阻力系数随雷诺数的降低而减小. 因此有必要研究微反应通道结构，通过优化流道内部结构，增强传热传质和提高反应效率. 受分形原理在微流控芯片产品方面成功应用的启发，对微型反应器管道结构尺寸优化，考虑到蛇形微反应器通道均匀分布，为了减少计算量，选择该微反应器的1/2作为计算域.

T型汇流微反应器中的蛇形通道内部流场分布比较复杂，是固相分布、入口条件、通道曲率变化等由多种因素共同作用的结果. 为了探究各种因素对该内部流场的演变机理，采用CFD-DEM耦合分析理论建立流体动力学模型，充分考虑颗粒受力情况（接触与非接触），能够准确地捕获颗粒运动轨迹及不同时刻的颗粒位置和速度.

固相颗粒采用DEM方法进行建模，固液两相间不存在质量传递，故建立流体相控制方程应考虑网格单元固体颗粒占据后流体相的体积分数 ${\alpha _{\rm{f}}}$、粒子间相互作用力 ${{{F}}_{\rm{a}}}$，计算流体单元对颗粒的平均作用力为

(1) ${{S}} = {{\sum\limits_i^n {{{{F}}_{\rm{a}}}} }}/{{{V_{\rm{f}}}}}.$

式中： ${{n}}$、 ${V_{\rm{f}}}$分别为单一流体网格单元内固体颗粒的数量和单个流体网格体积.

基于上述假设，流体控制方程修改为

(2) $\frac{\partial }{{\partial t}}\left( {{\alpha _{\rm{f}}}{\rho _{\rm{f}}}} \right) + \nabla \cdot \left( {{\alpha _{\rm{f}}}{\rho _{\rm{f}}}{{{\nu}} _{\rm{f}}}} \right) = 0,$

(3) $ \dfrac{\partial }{{\partial t}}\left( {{\alpha _{\rm{f}}}{\rho _{\rm{f}}}{{{\nu}} _{\rm{f}}}} \right) + \nabla \cdot \left( {{\alpha _{\rm{f}}}{\rho _{\rm{f}}}{{{\nu}} _{\rm{f}}}{{{\nu}} _{\rm{f}}}} \right) = - {\alpha _{\rm{f}}}\nabla p + \nabla \cdot { {{{{\tau}} _{\rm{f}}}}} + {\alpha _{\rm{f}}}{\rho _{\rm{f}}}{{g}} - {{S}}. $

式中： ${\rho _{{{\rm{f}}}}}$和 ${{{\nu}} _{{{\rm{f}}}}}$分别为流体相密度、速度； ${{{{{\tau}} _{{{\rm{f}}}}}}} $为液体相应力应变张量，

(4) $\begin{array}{l} { {{{{\tau}} _{\rm{f}}}}} = {\alpha _{\rm{f}}}{\eta _{\rm{f}}}\left( {\nabla {{{\nu}} _{\rm{f}}} + \nabla {{{\nu}} _{\rm{f}}}^{\rm{T}}} \right) + {\alpha _{\rm{f}}}\left( {\lambda {}_{\rm{f}} - \dfrac{2}{3}{\eta _{\rm{f}}}} \right)\nabla \cdot {{{\nu}} _{\rm{f}}}{I} . \end{array}$

其中 ${\eta _{\rm{f}}}$和 ${\lambda _{ {\rm{f}}}}$分别为流体相剪切黏度和体积黏度， $I $为流体相湍动能强度^[25].

通过计算上述流体相控制方程，能够得到微反应器流场中液相速度场及压力场，所求流场相解为固体颗粒运动提供了初始条件. 在计算CFD-DEM耦合的求解过程中，选择将固体颗粒从传统欧拉多相流方程中分离出来，利用DEM离散单元进行数值求解，在DEM框架下颗粒以拉格朗日方式追踪，动力学方程为

(5) ${m_i}\frac{{{\rm{d}}{{\nu}} {}_i}}{{{\rm{d}}t}} = {{{F}}_{\rm{a}}} + {{{F}}_{\rm{c}}} + {{G}},$

(6) ${I_i}\frac{{{\rm{d}}{{{\omega}} _i}}}{{{\rm{d}}t}} = {{{T}}_{\rm{c}}} .$

式中： ${{{F}}_{\rm{c}}}$和 ${{G}}$分别为颗粒碰撞接触力及颗粒本身重力. 颗粒在流场中的运动是受到多个作用力的结果，其中包括曳力 ${{{F}}_{\rm{d}}}$、压力梯度 ${{{F}}_{\rm{p}}}$、虚拟质量力 ${{{F}}_{\rm{v}} }$、Basset力 ${{{F}}_{\rm{b}}}$、Saffman力 ${{{F}}_{\rm{s}}}$、Magnus升力 ${{{F}}_{\rm{m}}}$及浮力 ${{{F}}_{\rm{f}}}$. 在化学反应进行过程中，液体密度远小于固体颗粒密度，可以忽略虚拟质量力及Basset力，可得

(7) ${{{F}}_{\rm{a}}} = {{{F}}_{\rm{d}}} + {{{F}}_{\rm{p}}} + {{{F}}_{\rm{s}}} + {{{F}}_{\rm{m}}} + {{{F}}_{\rm{f}}} .$

曳力求解方法为

(8) ${{{F}}_{\rm{d}}} = \frac{{{C_{\rm{d}}}}}{8}{\text{π}}{d_{\rm{p}}}^2{\rho _{\rm{p}}}\left| {{{{\nu}} _{\rm{f}}} - {{{\nu}} _{\rm{p}}}} \right|\left( {{{{\nu}} _{\rm{f}}} - {{{\nu}} _{\rm{p}}}} \right){\alpha _{\rm{f}}}^{ - (\chi + 1)},$

(9) $\chi = {\rm{3}}{\rm{.7 - 0}}{\rm{.65exp}}\;\left[ { - {{\left( {1.5 - \lg R{e_{\rm{p}}}} \right)}^2}/2} \right].$

式中： ${\rho _{\rm{p}}}$、 ${{{\nu}} _{\rm{p}}}$及 ${d_{\rm{p}}}$分别为固体颗粒密度、速度及直径； ${C_{\rm{d}}}$和 $R{e_{\rm{p}}}$分别为曳力系数以及颗粒雷诺数，

(10) ${C_{\rm{d}}} = {\left( {0.63 + \frac{{4.8}}{{R{e_{\rm{p}}}^{0.5}}}} \right)^2},$

(11) $R{e_{\rm{p}}} = \frac{{{\alpha _{\rm{f}}}{\rho _{\rm{p}}}{d_{\rm{p}}}\left| {{{{\nu}} _{\rm{f}}} - {{{\nu}} _{\rm{p}}}} \right|}}{{{\mu _{\rm{f}}}}},$

其中 ${\mu _{\rm{f}}}$为液体相动力黏度. 压力梯度、Saffman升力和Magnus升力分别由下式计算求得：

(12) ${{{F}}_{\rm{p}}} = \frac{1}{6}{\text{π}} {d_{\rm{p}}}^3\frac{{{\rm{d}}{{p}}}}{{{\rm{d}}t}},$

(13) ${{{F}}_{\rm{s}}} = 1.615d_{\rm{p}}^2\sqrt {{\rho _{\rm{f}}}{\mu _{\rm{f}}}} \left( {{{{\nu}} _{\rm{f}}} - {{{\nu}} _{\rm{p}}}} \right)\sqrt {\left| {\frac{{{\rm{d}}{{{\nu}} _{\rm{f}}}}}{{{\rm{d}}y}}} \right|},$

(14) ${{{F}}_{\rm{m}}} = \frac{\text{π}}{8}d_{\rm{p}}^3{\rho _{\rm{f}}}\left( {\omega \left( {{{{\nu}} _{\rm{p}}} - {{{\nu}} _{\rm{f}}}} \right)} \right)\left[ {1 + \theta \left( R \right)} \right].$

式中： $\theta \left( {{R}} \right)$为数量级较小的余量.

鉴于微反应器结构多管道及不同曲率半径等特点，管道内部流场动量方程应该添加雷诺应力梯度项. 为了保证方程封闭性求解，选择应用最广泛的可实现（Realizable） ${{k - }}\varepsilon $模型，运输方程为

(15) $\begin{split} \frac{\partial }{{\partial t}}\left( {\rho k} \right) + \frac{\partial }{{\partial {x_j}}}\left( {\rho {{\bar u}_j}k} \right) = & \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _{\rm{t}}}}}{{{\sigma _k}}}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right] + \\ & {\mu _{\rm{t}}}\left( {\frac{{\partial {{\bar u}_i}}}{{\partial {x_k}}} + \frac{{\partial {{\bar u}_{{k}}}}}{{\partial {x_i}}}} \right)\frac{{\partial {u_i}}}{{\partial {x_{{k}}}}} - \rho \varepsilon . \end{split}$

(16) $\begin{split} \frac{\partial }{{\partial t}}\left( {\rho \varepsilon } \right) + \frac{\partial }{{\partial {x_j}}}\left( {\rho {{\bar u}_j}k} \right) = & \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _{\rm{t}}}}}{{{\sigma _{\rm{\varepsilon }}}}}} \right)\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right] + \\ & {C_1}\rho E\varepsilon - {C_2}\rho \frac{{{\varepsilon ^2}}}{{k + \sqrt {\upsilon \varepsilon } }}. \end{split}$

式中：普朗特数 ${\sigma _k} = 1.0$， ${\sigma _{\rm{\varepsilon }}} = 1.2$，经验常数 ${C_{\rm{2}}} = 1.9$，

(17) ${C_{\rm{1}}} = {\rm{max}}\left\{ {0.43,\frac{\eta }{{\eta + 5}}} \right\}.$

其中 $\eta \!=\! {\left( {2{E_{ij}} \bullet {E_{ij}}} \right)^{1/2}({k}/{\varepsilon })}$， ${E_{ij}} \!=\! \left( {{{\partial {u_i}}}/{{\partial {u_j}}} + {{\partial {u_j}}}/{{\partial {u_i}}}} \right)/2$，湍动能黏度 $\,{\mu _{\rm{1}}} = \rho {C_{\rm{\mu }}}{{{k^2}}}/{\varepsilon }$，改进后的湍动黏度系数为

(18) ${C_\mu } = \frac{1}{{{A_0} + {A_{\rm{s}}}{U}k/\varepsilon }},$

${{{C}}_\mu }$与应变率有关，不是常数^[26].

考虑微反应器流道中的固体颗粒相互作用，采用Hertz-Mindlin无滑移接触软球模型进行计算；流动过程颗粒间会发生接触碰撞，颗粒碰撞受力方程为

(19) ${{{F}}_{\rm{c}}} = {{{F}}_{\rm{n}}} + {{{F}}_{\rm{t}}}.$

式中： ${{{F}}_{\rm{n}}}$为法向接触力，软球模型将颗粒在法向方向的接触过程用弹簧和阻尼器代替； ${{{F}}_{\rm{t}}}$为切向接触力，用弹簧、阻尼器和滑动器代替. 法向接触力方程、切向力方程分别为

(20) $\begin{array}{l} {F_{\rm{n}}} = - {K_{\rm{n}}}{d_{\rm{n}}} - {N_{\rm{n}}}{\nu _{\rm{n}}}. \end{array}$

(21) $\begin{array}{l} {F_{\rm{t}}} = \left\{ \begin{array}{l} - {K_{\rm{n}}}{d_{\rm{n}}} - {N_{\rm{t}}}{\nu _{\rm{t}}},\,\, \left| {{K_{\rm{t}}}{d_{\rm{t}}}} \right| <\left| {{K_{\rm{n}}}{d_{\rm{n}}}} \right|{C_{{\rm{fs}}}}; \\ \dfrac{{\left| {{K_{\rm{n}}}{d_{\rm{n}}}} \right|{C_{{\rm{fs}}}}{d_{\rm{t}}}}}{{\left| {{d_{\rm{t}}}} \right|}},\,\, \left| {{K_{\rm{t}}}{d_{\rm{t}}}} \right| \geqslant \left| {{K_{\rm{n}}}{d_{\rm{n}}}} \right|{C_{{\rm{fs}}}}. \end{array} \right. \end{array}$

弹簧标识的法向刚度 ${K_{\rm{n}}} = ({4}/{3}){E_{{\rm{eq}}}}\sqrt {{r_{{\rm{eq}}}}{d_{\rm{r}}}} $，阻尼器的法向阻尼 ${N_{\rm{n}}} = \sqrt {(5{K_{\rm{n}}}{m_{{\rm{eq}}}})} {N_{{\rm{ndomp}}}}$，切向阻尼 $ {N_{\rm{t}}} =$ $ \sqrt {(5{K_{\rm{t}}}{m_{{\rm{eq}}}})} {N_{{\rm{ndomp}}}}$，切向刚度 ${K_{\rm{t}}} = 8{G_{{\rm{eq}}}}\sqrt {{d_{\rm{n}}}{r_{{\rm{eq}}}}} $.

超声波是一种波动形式和能量形式的结合体，能够使流体获得动能加快流动，在一定程度上增强流场的湍动能，造成流场形成局部湍流，同时可能伴随空化产生，导致近壁面处局部高温高压，产生高速微射流和射击流等^[27]. 这种微气泡会因声场作用振荡、变大甚至崩溃的过程称为空化现象，这种空化效应促进反应的快速进行，从而提高产出率.

T型汇流下的蛇形微反应器管道内部是由颗粒与流体混合而成的黏性液固两相流，属于液固两相流范畴. 上述内容已论证液固两相流耦合相互作用，超声波技术的高频振动特性能够使流体中的流态参数发生变化，其周期性变化对流场调控造成影响，超声空化和湍流之间是相互激励、相互影响的. 在加载超声波激振的微反应器流场中，声场与固液两相流相耦合，会使流道内部的流场环境变得异常复杂，导致各种因素的物理分布情况很难直观分析. 基于对固液两相流数学模型的分析研究，结合超声波激振技术在流场产生的压力波对颗粒运动的影响.

在微型反应器加载超声波激振能够实现扰流，增强流场变化，达到抑制颗粒团聚现象. 在仿真模拟过程中，探头声压呈正弦变化，选择超声波频率为20 kHz，声压振幅为5；在T型流道交叉口Z轴垂直方向加入超声源，函数表达式为 $y = {p_{\rm{a}}}\sin\;(2{\text{π}} ft)$(其中p_a为声压振幅)；为了便于观察颗粒在流道内部的流动变化，选用非定常求解方式计算；超声波频率为20 kHz，周期为0.000 05 s；为了得到不同周期的参数值，将时间步长设为0.000 05 s，时间步数设为1 000，流体的物理运动时间为0. 005 s，足以达到稳定状态.

耦合超声波的T型汇流微反应器流道的内部流场分布比较复杂，颗粒群呈混沌态分布. 分形学试图透过混乱现象和不规则构型，揭示局部与整体的本质联系和特定规律. 在实际的生产或应用过程中，分形常与混沌紧密联系，混沌状态表示的是一种无规则性或者无周期性规律，混沌运动具有高度无序、混乱性、非线性等特点. 本文在探究耦合超声波技术对微反应器流场调控与抑制团聚现象的过程中，微反应器管道内部颗粒在流场中的运动轨迹近似一种混沌状态，很难实时捕捉. 本文描述颗粒运动轨迹的混沌状态并利用分形盒子算法研究流场，建立耦合超声波作用的动力学模型，模拟颗粒在管道内部的运动行为以及流场的演变机理，分析颗粒运动轨迹的混沌特性. 微反应器内部流场的盒维数计算方法为

(22) $\lg \left( {N\left( k \right)} \right) = \lg K + D\lg \left( {1/k} \right).$

式中：K为常数，k表示盒子的尺寸，N（k）为尺寸为k条件下覆盖所有图像像素的最少盒子数.

面向微反应器流场识别的分形处理流程如图2所示. 在流场图片获取后，先进行灰度转换、特征检测、边缘提取等操作，提取目标区域；然后采用分割阈值方法进行二值化处理，转换为黑白位图；再利用数学形态方法去除噪声，最后将图片放入到盒维数算法程序中，实现分形维数的计算.

Pearson等^[28]的研究表明，湍流的旋涡嵌套结构具有分形自相似特征，不同流型的分维数之间存在差异，可以作为流型区分的依据. 本文中微反应器的流场流动、颗粒运动是随机可变的，具有显著非线性特征. 此外，颗粒运动实质上是一种非线性耗散的运动系统，表现出一种混沌状态，具有分形性质. 超声波传输机理与颗粒的相互作用影响是非常复杂的，基于分形盒子算法研究超声波激振对管道内部颗粒集群的分形特征，对比分析有、无超声波加载微反应器交汇处以及管道处的分形维数. 由于超声波具有分行标度不变性，计盒维数能够反映超声激振作用后微反应器内部颗粒集群的运动轨迹和分布特点.

基于提出的CFD-DEM建模方法，可以建立微反应器流场动力学模型. 微反应器采用圆形截面通道，直径为1 mm，选择液态水和SiC颗粒分别作为液相和固相进行模拟，其中颗粒总数目为40 000，直径为40 μm，密度为2 500 kg/m³，泊松比、恢复系数、静摩擦系数、动摩擦系数分别为0.25、0.90、0.15、0.1.

微反应器蛇形管道流场分布通过CFD方法获取，选择不可压缩的流体、可实现 $k{\rm{ - }}\varepsilon $紊流模型求解；对于固相颗粒，在DEM计算中采用拉格朗日坐标方法，分析颗粒及颗粒集群在管道内部的运动特性和分布规律. 基于上述假设可知，基于CFD-DEM的微反应器流场数值计算流程如图3所示. 首先通过CFD模拟流场流动特性，然后将流场数据传输到DEM模型中；DEM计算流场中的颗粒受力及颗粒集群流动轨迹，通过内部迭代计算颗粒下一时间步长运动并实时更新位置和速度，最后将数据结果再次传输到CFD中，进行下一时间步长迭代计算.

为了获得较好的收敛性和计算准确性，考虑到双入口通道的分叉结构及边界层流，为了保证计算精度对汇流交汇处结构局部网格加密，计算结果的准确性不受网格密度的影响，采用四面体和六面体网格进行划分模拟. 数值计算与实验值之间的误差来源主要有：物理模型近似误差、差分方程阶段误差、求解区域的离散/迭代/舍入误差等（可以参考专家一意见第一点）. 在计算过程中，离散误差随网格变细而变小，但由于网格变细时，离散点数增多，舍入误差随之加大.在工程流体计算中，不规则区域内如何生成网格对计算结果及效率具有重要的影响. 本文建模时的网格划分采用结构化网格并局部加密，网格质量在计算条件允许的范围内调整到最高，以保证计算结果的精确性与可重复性；在此基础上，完成了网格无关性验证，通过图4可以看出，3 种网格划分方式在相同条件下得到的计算结果具有很高的重合度，误差约为1.6%，确保本文的数值计算求解均为网格独立解.

对CFD-DEM模型边界设置，该模型采用有限体积法对边界条件作用下的控制方程进行离散. 通道双入口边界设置为速度入口，分别设定流体与固体颗粒的初始速度以及颗粒体积分数，保证进入管道内部前颗粒与流体进行充分混合. 考虑到颗粒尺度小，固体颗粒与流体速度须根据情况来调整.

通道出口边界设置为压力出口，选择出口处作为压强参考点，数值为一个标准大气压. 此外，微反应器液-固交界的壁面设定为速度无滑移表面，考虑管道轴向壁面对颗粒的作用力. T型管道交汇处垂直方向引入超声场，实现声场与流场、声场与固相等多物理场耦合. 采用SIMPLEC算法实现该微反应器内部压力与速度的耦合，非稳态项采用二阶隐式，动量方程采用二阶迎风离散，调整各参数松弛因子，缩短收敛时间.

在数值模拟过程中，当流道内流场分布达到稳定后，向通道内加入反应动力学模型，模拟通道内部流动碰撞、化学反应过程. 为了确保计算的准确性，当各参数残差小于10⁻⁶且能量残差小于10⁻¹²时进出口质量守恒和能量守恒，则认为模拟结果收敛.

基于多组不同工况参数的CFD-DEM数值模拟结果，对有、无超声波作用下速度、湍动能、压力等参数进行对比分析. 结合可视化软件监测颗粒的运动轨迹，得到小尺度颗粒接触碰撞的运动规律，归纳出微反应器内部流场演变机理和混合强度，为T型汇流微反应器的流场演化和破团聚研究提供理论参考.

考虑超声波激振对流场分布的影响，选取4组不同速度，作为研究T型汇流微反应器蛇形流道流场演化的入口条件. 微反应器的2个速度入口分别是沿x轴和沿y轴方向，z轴表示管道的竖直高度方向. 因多相流数学模型剖面图能够较完善反映不同速度作用下流场的变化规律，选取几何模型截面观测微反应器管道内部流动情况及颗粒集群团聚位置，揭示流体流动演变机理和受扰动内部流场的变化规律. 利用超声波激振技术产生波动能量来驱动和扰动流体，抑制甚至破除团聚现象，实现流体的均匀离散分布. 加载超声波前、后的速度场、动压力曲线如图5所示.

如图5（a）所示为加载超声波作用下的速度运动特性曲线. 4组不同速度入口工况下的流场截面速度变化均呈正弦波动，符合能量波动规律. 速度为2 m/s时波峰和波谷保持在6和2左右；2.5 m/s时波峰和波谷保持在7和3左右；3 m/s时波峰和波谷分为2段现象：0~1.0 s时波峰和波谷维持在9和4.5左右，1.0~2.0 s时波峰和波谷维持在11和6左右；3.5 m/s时波峰和波谷维持在9和5左右. 由图5可知，当速度为3 m/s流动时，微反应器管道截面速度峰值较大，说明在同种超声波频率作用下，微反应器管道入口速度在3 m/s时受超声波动效果显著有利于加速溶液混合流动，速度压力曲线均能体现出. 3 m/s与3.5 m/s二者速度和动压变化趋势顺序不同是受到波激振影响的缘故. 对比有、无超声场可以看出，超声波能够改变原速度场分布，在一定范围内增大速度，增强湍流强度；否则对于微型管道而言会引起交汇处溶液混合部分颗粒获取较大速度，朝着不同方向运动，出现颗粒回流现象. 图5（b）未加载超声波激振，对比前者4组速度出现不同程度的抖动，仅能够证明无超声波激振作用下的微反应器管道内部流场流动为层流，该层流流动现象与文献[29]的现象相似，无法准确获取管道内部的流场流动规律；由于流速波动性弱，可能导致管道弯折处出现颗粒团聚现象，在0.3 s后速度波动趋于平缓. 对比以上结果可以发现，在加载超声波作用下的微反应器速度波动性大，增强流体流动有利于加快混合. 如图5（c）、4（d）所示分别为加载超声波作用和未加载超声波作用的截面的动压变化曲线. 针对加载超声波作用下的管道截面动压呈超声波激励正弦波动，在一定范围内随着速度的增大，动压波动增大. 对比3 m/s和3.5 m/s两组速度可以看出，在0.7 s后，前者速度峰值略高于后者，说明3 m/s的速度条件适合进行混合.

1）流场特征参数云图. 鉴于上述的流场速度、动压变化曲线，仅能够从整体上观察微反应器内部截面的变化趋势. 为了更清楚地认识流场的流动情况，掌握更多的流场细节数据，选取不同入口速度下某一时刻的速度、动压云图，如图5、6所示.

由图5的速度、动压变化曲线可知，流场变化在0.4 s左右趋于平缓. 为了更好地观察流场流动情况，选择1 s时刻进行监测. 根据图6可以看出，4种工况参数下流场的速度变化差异不是太大. 在流道交汇处出现融合现象，由于管道尺寸与曲率半径改变，此处会造成催化剂颗粒碰撞加剧，能量损失，局部有一定的漩涡形成汇流. 在交汇处后与管道弯折处之间的部分，混合溶液流动平缓，说明该段范围内混合流动弱. 管道弯折处由于流动方向改变和曲率半径减小，携带固体颗粒的溶液在此处会造成管道阻塞，从管道弯折处内部与外部出现速度云图. 不规则的云图表示速度在此处的改变. 根据图7的动压云图变化可以看出，流场混合溶液压力在T型管道入口处较小，经过交汇处内部动压力逐渐增大，管道边缘处动压小于流场管道内部，在管道弯折处由外到内动压力呈现梯度变化，这是由于颗粒的汇聚导致内管道阻塞，混合溶液流通不畅，不考虑液体的黏性，忽略黏性阻力因素. 针对该微反应器管道阻塞、反应速率缓慢等问题，引入超声波激振技术对流场速度、动压力等工况参数进行流场调控，速度与动压力分布如图8、9所示.

根据图8、9可知，速度云图、动压力云图在超声波作用下对流场调控均匀分布，说明超声波激振能够优化管道内部流场分布. 从图7可以看出，4种速度条件在超声波激振作用下，截面速度云图分布均匀，弯折处未出现急剧变化区域，说明此处溶液流动均匀分布无颗粒汇聚情况. 通过图8所示的结果可以发现，动压在超声波激振作用下管道整体分布均匀，无明显梯度变化，说明超声波对于动压场分布具有一定的调控作用.

2）多组超声波频率下流场演变规律.

图6~9所示结果的加载超声波激振频率为20 kHz，未涉及更高频率超声波对流场的影响和演变规律. 为了考察不同频率激振对微反应器流场参数的影响，提供了4组激振频率条件下的流场分布结果，如图10所示.

由图10（a）可以看出，超声波频率取25 kHz，速度曲线波动幅值较接近；超声波频率为30 kHz时的速度正弦波动稍弱；15 kHz作用下在1.5 s后速度波动较三者有明显的变化；从峰值和能量消耗及成本因素考虑，频率为20 kHz在4组频率中性能较稳. 如图9（b）所示为4组不同频率作用下的动压力变化. 通过对比4组波动曲线可以看出，30 kHz动压力曲线波动较弱，相对于其他3组频率变化效果不明显，说明该条件下动压变化不平缓，流场内部流体混合不剧烈，化学反应进程缓慢. 当超声波频率为25、15 kHz时，动压力、速度非均匀变化不利于流场颗粒充分接触，不利于催化反应进行. 由于频率越高，成本花费较大且得到的效果不显著，综合考虑选择超声波频率为20 kHz较合适.

颗粒碰撞是一种接近混沌状态的流固耦合，研究颗粒-颗粒、颗粒-壁面接触碰撞情况，但该系统具有高度非线性、不确定性及复杂性，目前对这方面的研究较少. 本文以混沌和分形方法为理论基础应用处理非线性、复杂性的多相流系统，研究微反应器内部的拓扑结构，揭示颗粒运动轨迹及接触碰撞规律.

未加载超声波条件下的微反应器管道颗粒分布如图11所示. 由于入口汇流作用，蛇形微反应通道内部流场流动分布不均匀，T型微反应器管道弯折处和交汇处颗粒汇集现象明显. 为了考察弯道处的颗粒分布细节，选择第一个弯道局部剖切，可以看出图11（a）~（d）都出现颗粒的汇聚现象，伴随速度的增大，颗粒集中汇聚在管道外侧，内侧及中间部分颗粒分布较少；入口速度越大，汇聚现象越明显. 从图11（a）~（d）可知，从管道颗粒的流动分布和流向情况能够较直观地看出颗粒集群团聚位置. 颗粒在微反应器流动进行催化反应会在管道处出现团聚现象，导致管道阻塞不利于该反应的进行. 由管径的变化可以看出，颗粒直径突变处存在颗粒集聚，颗粒之间明显形成团聚，影响管道内部流体的运动，严重影响混合反应，降低反应效率.

加载超声波条件下的微反应器管道颗粒分布如图12所示. 在加载超声波激振作用下，微反应器内部流动分布均匀，随着速度的增加管道处颗粒团聚现象得到抑制，可以验证加载超声波有利于抑制颗粒团聚，有利于催化反应的进行. 从图12（a）~（d）可知，加载超声波后颗粒分布更加均匀，弯道颗粒流动较稳定；对比无超声波管道剖面颗粒分布可以看出，管道内部流动平缓且分布均匀，对混合溶液和流场的调控具有积极意义.

1）颗粒集群参数的分布特征.

在固液两相流动过程中，相界面分布和形状是随着流场流动变化的，多相流流动过程较单相流复杂，比单相流多一些新特性：流型、相含率、压降等. 通过对截面的含固率统计测量，可以得出两相流中分散相的浓度与分布规律. 由于无超声波作用的微反应器内部管道流动或催化反应过程中出现颗粒团聚问题，会导致管道阻塞，影响化学反应速率，降低产出率等. 对于管道内部的阻塞状况，难以进行直接观测，需要借助数值计算间接手段得到流场特征.

针对颗粒集群的运动特征，选取截面上19个样本点进行监测，结果如图13、14所示.

对比图13（a）、（b）可以看出，未加载超声波作用下颗粒在3个弯折峰值达到0.35以上，说明该位置出现颗粒团聚现象；在加载超声波作用后，样本点颗粒体积分数下降到0.35以下，而其他位置处颗粒分布基本无明显变化，能够说明加载超声波后颗粒在弯道处分布密度降低，抑制颗粒阻塞管道造成团聚；超声波激振有利于微反应器管道内部颗粒分布，也有利于流场的分布与调控；在入口速度为2 m/s的条件下，在第7、17特征点的颗粒体积分数相比未加载超声波作用下降低；在入口速度为3.5 m/s的条件下，颗粒体积分数下降更加明显.

前面涉及到特征点位置的颗粒体积分数分布，特征点位置速度曲线与颗粒集群体积分数具有相关性. 根据图14的特征点位置速度分布曲线可以发现，3个弯道处颗粒体积分数汇聚达到30%~40%，相反速度在这些位置减弱到1~3 m/s. 弯道处颗粒大量汇聚后，在一定程度上导致流速减小，进而影响流场流动. 在图14（b）的加载超声波作用后，特征点位置速度均增大，说明超声波激振能够在一定程度上调控流场流速和优化流道分布. 随着入口速度的增大，特征点位置流速和湍动能增大，有利于抑制颗粒团聚现象.

2）颗粒集群的分形特征.

考虑微反应器流道中颗粒集群分布的分形特征提取问题，可以根据不同工况下的颗粒运动变化，将获取到的图片数字化处理后求解分形维数. 通过仿真T型汇流下蛇形微反应器管道内部催化剂的颗粒分布，选取有无超声波2组T型交汇处、第一个管道位置流道，开展研究颗粒集分布的计盒维数计算.

通过对未加载超声波激振的微反应器内部管道进行监测，针对交汇处、弯道处颗粒集群进行捕捉，利用图像处理技术进行图像预处理；然后利用盒维数方法计算不同参数条件作用下的分形维数D，结果如图15、16所示. 不同速度条件作用下的微反应器内部交汇处、弯道处颗粒集群的分形维数分布不同，由此可以得出作用条件不同，分形维数不同，微观形貌存在差异，但基本都保持在2.0左右. 实际上，图像的计盒维数越接近2，表明所关心的区域越能够反映整体的图像范围分布；计盒维数越接近1，表明所关心的区域与整体区域趋于线性分布. 由图15、16可知每组条件作用下的颗粒集群分布，分形维数分别为2.15和2.13、1.98和2.12、2.02和2.13、2.04和2.12，表明此时流道内部颗粒集群具有分形特征，已经达到混沌状态，也说明交汇处和弯道处已达到混沌状态.

加载超声波条件下的交汇处与弯道处的颗粒分布分形维数如图17、18所示. 微反应器管道交汇处和弯折处的颗粒分布分形维数分别是1.95和1.81、2.08和1.82、2.03和2.05、2.07和2.06. 通过对比图15、16可以看出，是否耦合超声场对分形盒维数的影响不大，验证了超声波保持分形标度的不变性.

（1）在同种超声波频率的作用下，微反应器管道进口速度在一定范围内增大，有利于加速混合溶液流动，提高反应效率.

（2）流道交汇处出现融合现象，由于管道尺寸与曲率半径改变，会造成固体颗粒碰撞加剧，能量损失.

（3）动压在超声波激振作用下分布均匀，局部区域出现空泡现象，在一定程度上影响湍流强度，对于提高混合反应效率具有正向作用.

（4）超声波激振能够在一定程度上调控流场流速和优化流道分布，随着入口速度的增大，特征点位置流速和湍动能增大，有利于抑制颗粒团聚现象.

（5）针对交汇处、弯道处的颗粒集群进行捕捉，利用盒维数方法计算不同参数条件作用下的分形维数. 结果表明，此时流道内部颗粒集群具有分形特征，处于混沌状态.

（6）耦合超声场对分形盒维数的影响不大，验证了超声波保持分形标度的不变性.

声场耦合条件下的微反应器多相流建模与流场特性分析具有较高难度，本文在该方面进行了有益尝试，可以为微反应器生产工艺优化及燃料电池、微流控芯片等领域的机理建模与应用系统研发提供理论依据与技术支持. 后续研究工作将围绕微反应器内部流场跨尺度建模与振动动力学特性展开.