填充床反应器中的有限时间混合与反应

《Journal of Fluid Mechanics》：Finite-time mixing and reaction in packed-bed reactors

　　摘要：我们研究了有限时间内多孔介质中的混合动力学，采用了结合拉格朗日粒子追踪的孔隙尺度格子-Boltzmann模拟方法。我们基于Lester等人（2018年《流体力学杂志》855卷第770-803页）提出的移动Protean框架方法，计算了随机堆积床中流体的变形。从提取的拉格朗日运动学数据中，我们构建了一个基于层状聚集的混合模型，该模型能够很好地预测模拟得到的欧拉标量场。我们的研究结果揭示了早期混合阶段的特征：这一阶段主要由剪切驱动的流体变形主导，溶质混合源于扩散浓度元素的随机重叠。在这个阶段，混合过程进展缓慢，并且浓度方差随时间衰减，遵循公式$\sigma _c^2 \propto \textit{Pe}^{-\alpha /(2\alpha +1)} t^{-1/2}$，其中$\textit{Pe}$是Péclet数，$\alpha$是表征剪切变形的指数。这种动态现象发生在Péclet数相对较小，且与指数混合时间尺度和剪切变形时间尺度之比较小时。分析还表明，在有限尺寸的随机堆积床中，剪切诱导的混合过程控制着流体-固体界面上早期反应的均匀化，这类系统的Péclet数通常在$ \textit{Pe} =O(10^2)$的范围内。

1. 引言：尽管多孔介质中的完全混合机制已经被广泛研究，但在早期有限时间内的混合动力学仍然了解不多。在尺寸有限或以有限时间过程（如反应）为特征的多孔系统中，有限时间内的混合在控制注入溶液的均匀化方面起着决定性作用。例如，在小型紧凑的多孔反应器中，尽管分子扩散速度较慢，不足以支持混合，但孔隙尺度上的对流有助于在反应器的纵向尺寸$L$范围内有效分布溶质。然而，停留时间$L/U$必须足够长，以确保分子的完全转化（参考文献Maggiolo等人，2020年）。这意味着设备必须在中等Péclet数$ \textit{Pe} =O(10^2)$的条件下运行。因此，对流驱动的混合效果较弱，并受到有限物理反应器空间$L$的限制。类似的小型多孔系统还应用于其他领域，如膜和电极（参考文献Farzaneh等人，2021年）。此外，在高反应性介质中，由于反应的抑制，混合可能在非常短的时间内就受到限制。在这些条件下，可能无法实现完全混合，混合动力学在有限时间内才能完全体现。无论是因为反应时间快还是介质物理尺寸短，溶质混合所需的时间可能都短于标志渐近行为开始的特征混合时间$t_m$。这个混合时间与特征孔隙长度$d$成正比，与平均流速$U$成反比，即$t_m \propto d/U$。更准确地说，$t_m$反映了由孔隙尺度流动和几何形状决定的特定混合机制。在Péclet数有限或较大时，通过多孔介质传输的溶质往往会形成由低浓度空隙分隔的细长高浓度细丝（参考文献Dentz和de Barros，2015年）。混合指的是标量场（如浓度）从最初的非均匀分离状态转变为最终的空间均匀分布（参考文献Villermaux，2019年）。在多孔介质中，这一过程依赖于流动在极早期通过孔隙尺度几何形状产生的运动学变形来消除小尺度浓度差异的能力。换句话说，混合完全由孔隙尺度结构介导。与二维流动不同，三维多孔结构本质上会产生混沌流动（参考文献Lester等人，2013年；Lester和Denz以及Le Borgne，2016年），其中随机发生的拉伸和折叠事件导致浓度场的长期指数级拉伸（参考文献Heyman等人，2020年）。无限小流体元素长度的平均长期增长率由Lyapunov指数$\lambda$量化，它决定了混沌混合的开始时间$t_m \sim d/(2 \lambda U) \ln (\textit{Pe})$。Turuban等人（2019年）进一步证明，在具有简单人工晶格结构的颗粒多孔介质中，尽管混沌混合的效果较弱，但它强烈依赖于平均流动方向相对于晶格几何形状的取向。对于随机拉伸取向的特定情况，Lyapunov指数被理论估计为$\lambda \approx 0.11$（参考文献Lester等人，2013年；Lester和Denz以及Le Borgne，2016年）。在随机堆积床中的实验研究报道了稍高的值：$\lambda = 0.20$（参考文献Kree和Villermaux，2017年），$\lambda = 0.21$（参考文献Heyman等人，2020年）和$\lambda = 0.17$（参考文献Heyman等人，2021年）。在反应时间较快或介质物理尺寸较短的情况下，溶质混合的时间可能短于标志渐近行为开始的特征混合时间$t_m$。这个混合时间与特征孔隙长度$d$成正比，与平均流速$U$成反比，即$t_m \propto d/U$。更确切地说，$t_m$反映了由孔隙尺度流动和几何形状决定的特定混合机制。在Péclet数有限或较大的情况下，通过多孔介质传输的溶质常常会组织成由低浓度空隙分隔的细长高浓度长丝（参考文献Dentz和de Barros，2015年）。混合指的是标量场（如浓度）从最初的非均匀分离状态转变为最终的空间均匀分布（参考文献Villermaux，2019年）。在多孔介质中，这一过程依赖于流动在非常早期通过孔隙尺度几何形状产生的运动学变形来消除小尺度浓度差异的能力。换句话说，混合完全由孔隙尺度结构介导。与二维流动不同，三维多孔结构本质上会产生混沌流动（参考文献Lester等人，2013年；Lester和Denz以及Le Borgne，2016年），其中随机发生的拉伸和折叠事件导致浓度场的长期指数级拉伸（参考文献Heyman等人，2020年）。无穷小流体元素长度的平均长期增长率由Lyapunov指数$\lambda$量化，它决定了混沌混合的开始时间$t_m \sim d/(2 \lambda U) \ln (\textit{Pe})$。Turuban等人（2019年）还证明，在具有简单人工晶格结构的颗粒多孔介质中，即使混沌混合效果较弱，它也强烈依赖于平均流动方向相对于晶格几何形状的取向。对于随机拉伸取向的特定情况，Lyapunov指数被理论估计为$\lambda \approx 0.11$（参考文献Lester等人，2013年；Lester和Denz以及Le Borgne，2016年）。在快速反应时间或多孔床物理长度仅限于几个颗粒直径的情况下，由于停留时间$1/k_r$或$L/U$（其中$k_r$分别为反应速率）与混合时间$d/(\lambda U)$相当，因此无法实现完全混合。Lyapunov指数是一个无限时间的度量，无法完全捕捉三维变形。在如此短的时间内，剪切诱导的变形和混沌混合共存。Lester等人（2018年）表明，在有限时间内，剪切诱导的变形可能占主导。Heyman等人（2021年）记录了在初始瞬态阶段之后混沌混合特征的显现，在此期间浓度差异仍然存在。这个瞬态阶段与混沌混合需要时间来完全发展是一致的（至少需要达到与混合的指数模式相对应的Batchelor尺度的时间）。Aquino等人（2023年）展示了剪切变形如何控制多孔介质中物种的长时间分布以及到达固体反应表面的旅行时间统计，从而有效地决定了多孔介质的反应效率。然而，导致这一渐近极限的瞬态阶段尚未得到充分探索。层状理论已被开发出来捕捉这种动力学。这些理论最初是为湍流提出的（参考文献Duplat和Villermaux，2008年），后来扩展到随机多孔流动（参考文献Dentz和Hidalgo以及Lester，2023年），这些模型将混合描述为高浓度层状结构（或片状结构）的形成及其随后的随机聚集。这个框架突出了一个早期的拉伸-压缩阶段，在此阶段，对流（i）在空间中随机重新分布高浓度区域，使先前分离的口袋靠近，并（ii）放大层状结构内的横向浓度梯度，从而加速扩散驱动的稀释。这两种过程都是强大的混合机制。在$t_m$之后，重叠的层状结构可以被建模为经历随机聚集，这种聚集可能是完全相关的，也可能是部分相关的，这取决于孔隙尺度的动力学（参考文献Duplat和Jouary以及Villermaux，2010年；Heyman等人，2024年）。因此，层状模型为研究有限时间内的混合提供了一个有用的框架，因为它也捕捉了早期流体诱导的变形动力学。这种动力学取决于几个因素，最显著的是（i）搅拌协议的空间维度和（ii）流动的分布和随机性。Villermaux（2019年）回顾了各种搅拌协议，包括二维剪切诱导的变形。对于这种情况，层状理论的预测准确性已经通过实验得到验证（参考文献Souzy等人，2018年），并且混合时间被证明与Péclet数有弱代数依赖性，即$t_m \sim (d/U) \textit{Pe}^{1/3}$。类似的缩放也适用于涡流诱导的混合（参考文献Villermaux，2019年）。在本工作中，我们利用计算分析在欧拉和拉格朗日框架内定量研究了单分散堆积床中的有限时间混合行为。我们解决的研究问题是：（i）在堆积床多孔介质中，哪种混合模式在有限时间内最为相关；（ii）这种混合模式如何影响流体-固体界面的孔隙尺度浓度分布；以及（iii）是否可以建模和预测有限尺寸堆积床反应器中表面反应的均匀性。从基本角度来看，解决这些问题对于深入理解多孔介质中的早期混合至关重要。从技术角度来看，需要新的模型来捕捉堆积床中流体-固体界面处物种的早期重新分布，从而改进小型多孔反应器的设计和操作控制。通常提出两种数值方法来量化孔隙尺度混合并提取相应的宏观反应速率。欧拉模型跟踪浓度场的时间演化（参考文献Gurung和Ginn，2020年），而拉格朗日方法，如随机反应行走（参考文献Valocchi、Bolster和Werth，2019年；Sole-Mari等人，2021年），则跟踪单个示踪剂的轨迹。这两种方法都面临计算挑战：欧拉模型需要高分辨率的孔隙尺度离散化，而拉格朗日模型受到可解析示踪剂数量的限制。其他建模框架，包括分形导数扩散模型（参考文献Liang等人，2023年）和有效反应传输方程（参考文献Dentz和Gouze以及Carrera，2011年），可以用来解释异质性、异常扩散和小尺度混合。然而，这些模型依赖于描述孔隙几何形状的微观结构参数，这些参数往往难以事先测量，使得直接的孔隙尺度模拟不可或缺。为了解决欧拉流动和孔隙尺度上的浓度场，我们使用了格子-Boltzmann方法（LBM），该方法非常适合解决复杂几何形状中的层流，因为它允许轻松地局部表示投影到规则空间网格上的孔隙表面。我们以有限长度的管状通道中包含的堆积床的几何形状作为测试案例。然后我们应用拉格朗日粒子追踪（LPT）算法来计算被流动携带的无质量颗粒（示踪剂）的拉格朗日轨迹。沿着这些轨迹，我们基于Lester等人的移动Protean框架（参考文献Lester、Dentz、Le Borgne和de Barros，2018年）计算变形张量，这最终使我们能够确定有限时间内流体诱导的主导变形模式。然后我们使用在计算出的变形上构建的层状模型来预测不同浓度区域的空间重组和混合。我们通过将层状模型的预测与欧拉浓度场的解进行比较来评估其有效性。本文的结构如下：第2节描述了数值方法，第3.1节和第3.2节研究了孔隙尺度速度和变形的特征，第3.3节和第3.4节分析了浓度场的演化，而第3.5节和第3.6节使用层状聚集模型解释了混合动力学，我们将在第3.8节将其扩展到反应性多孔系统。第4节总结了主要发现。

2. 数值方法
2.1. 几何构造
为了研究有限时间内的混合，我们以一个小的固定堆积床为例：准球形催化颗粒紧密堆积在管状反应器中（图1）。不断注入分子溶液，同时催化剂颗粒表面触发化学反应，逐渐转化溶质。一个显著的例子是分子太阳能热能存储（MOST）系统，其中光激活分子以化学键的形式储存太阳能（参考文献Wang等人，2019年；Maggiolo等人，2024年）。溶质浓度通常保持较低，以防止分子抑制，而堆积床的多孔几何形状为反应提供了高比表面积。选择小型反应器有以下几个优势：（i）能量提取效率高，（ii）催化剂的使用和维护成本低，以及（iii）可以紧凑地集成到小容器中。这些系统对孔隙尺度上的混合演化非常敏感。实际上，假设混合完全发展并且混合条件迅速发展可能会导致预测反应器性能时出现显著误差（参考文献Sole-Mari、Bolster和Fernàndez-Garcia，2023年）。图1。（a）两个半径为$R$的管状反应器，其中包含由直径为$d= R/3.1$（上图）和$d=R/4.7$（下图）的催化剂球形颗粒组成的反应性多孔介质。这两个系统被用作格子玻尔兹曼（lattice-Boltzmann）模拟的几何输入。(b) 标准化流速${u_{\kern-1pt j}^*(\boldsymbol{x})=u_{\kern-1pt j}(\boldsymbol{x})/U$（其中$U$为平均流向流速）的概率分布函数，沿着流向$u_x$（黑点标记）和横向$u_t=(u_y^2+u_z^2)^{1/2}$（虚线）方向；颗粒大小分别为$d= R/3.1$（上图）和$d=R/4.7$（下图）。我们使用开源计算机图形软件工具Blender来构建反应器几何形状。通过模拟球形颗粒落入圆柱形容器中的刚体动力学，来模仿由球形颗粒组成的催化剂反应器的微结构。我们考虑了两种几何配置，反应器半径与颗粒直径的比例分别为$R/d=3.1$和$4.7$。反应器长度相应地进行缩放，分别为$L=28d$和$42d$。颗粒沉积后，我们提取颗粒中心位置，并重建所得几何形状的二值化三维图像矩阵。这些二值化图像包含$7.5 \ 10^7$和$13.1 \ 10^8$个体素（voxels），体素与颗粒大小的比分别为$\Delta x/d \approx 1/40$和$\Delta x/d \approx 1/32$，对应$R/d=3.1$和$4.7$的配置。这种分辨率确保了用于后续数值格子玻尔兹曼计算的体素大小，始终小于与目标流速相关的最小特征传输尺度（在线性剪切下的Batchelor长度尺度$s_b/d \sim {\textit{Pe}}^{-1/3} \approx 1/7$，在填充分布床中的指数拉伸下$s_b/d \sim (0.1 \textit{Pe})^{-1/2} \approx 1/5$，详见§ 3.6）。这两种配置如图1(a)所示。

2.2. 孔隙尺度模拟
格子玻尔兹曼（LBM）被应用于解决孔隙尺度上的动量方程和溶质反应问题。由于其算法并行化的便利性，LBM特别适合这项任务，能够在高性能计算基础设施上实现高效的数值求解，包括复杂几何形状中的流动模拟。我们使用开源代码lbdm（Maggiolo 2024），该代码之前已被用于模拟多孔介质中的溶质传输和反应。Maggiolo、Modin和Kalagasidis（2023）提供了该算法在多孔介质应用中的验证。首先在多孔几何结构内求解动量方程，通过计算由模拟压力梯度的体力驱动的格子群体的流动和碰撞来实现。在催化剂截面的入口和出口添加了长度为$R$的两个缓冲区。这些区域是没有障碍物的反应器容器的延伸部分，允许流动沿流向重新对齐，并实现周期性边界条件。这种策略允许应用模拟压力梯度的体力，避免了边界压力的调整，提高了数值稳定性（Yang和Wang 2023）。在流体-固体边界处施加无滑移条件。达到稳态后，所有模拟中的流动Péclet数为$\textit{Pe}=Ud/D_m\approx 280$，其中$U$为平均流速，$D_m$为分子扩散率。在这种$\textit{Pe}$值下，预期会出现混沌混合特征（Lester、Metcalfe和Rudman 2014）。通过单松弛时间格子玻尔兹曼方程（2.1）求解流场：
$$
f_r (\boldsymbol{x} + \boldsymbol{c}_r, t+1) - f_r(\boldsymbol{x},t) = - \tau _\nu ^{-1} \big( f_r(\boldsymbol{x},t)-f^{eq}_r (\boldsymbol{x},t) \big) + F_r,
$$
其中$f_r(\boldsymbol{x},t)$表示位置$\boldsymbol{x}=(x,y,z)$和时间$t$沿$r$方向的分布函数，$\boldsymbol{c}_r$表示沿$r$方向的离散速度向量，$\tau _\nu$是松弛时间（与流体粘度成正比）。用$f_r^{eq}$表示沿$r$方向的平衡分布函数（2.2）：
$$
f^{eq}_r (\boldsymbol x,t) = w_r \rho \ \Bigg ( 1 +\frac {c_{\textit{rj}} u_{\kern-1pt j}(\boldsymbol x,t)}{c_s^2} + \frac {(c_{\textit{rj}} u_{\kern-1pt j}(\boldsymbol x,t))^2}{c_s^4} - \frac {u_{\kern-1pt j}^2(\boldsymbol x,t)}{2c_s^2} \Bigg ),
$$
其中$c_s$是声速，$w_r$是沿$r$方向的三维格子结构的19速格子权重参数。通过在每个计算单元积分分布函数的流体动力矩来计算流体场。然后得到稳态速度向量$u_{\kern-1pt j}(\boldsymbol{x})$和密度$\rho (\boldsymbol{x})$（2.3）：
$$
\rho (\boldsymbol x) = \sum _r f_r(\boldsymbol x), \nonumber
$$
$$
\rho (\boldsymbol x) {u}_j (\boldsymbol x) = \sum _r {c}_{\textit{rj}} f_r(\boldsymbol x) + \frac {1}{2}\bigg ( \frac {\Delta P}{L} \bigg )_{\kern-1.5pt j}.
$$
这里，${\Delta P/}{L}$是驱动流体通过介质的施加压力梯度。这一项也出现在体力$F_r$中（2.1）（Guo、Zheng和Shi 2002）：
$$
F_r (\mathbf{x},t)= \Bigg ( 1-\frac {1}{2\tau _\nu } \Bigg ) w_r \Bigg ( \frac {{c}_{\textit{rj}}-{u}_j(\mathbf{x},t)}{c_s^2} +\frac {{c}_{\textit{rj}}{u}_j(\mathbf{x},t)}{c_s^4}{c}_{\textit{rj}} \Bigg ) \ \bigg ( \frac {\Delta P}{L} \bigg )_{\kern-1.5pt r}.
$$
下一步，使用计算出的稳态流场通过第二个格子群体来解决被动标量浓度$c$的传输问题。标量在入口处以浓度$c_0$连续注入，并可能在催化剂-流体界面上发生反应。系统受稳态不可压缩动量方程（2.6）和控制：
$$
u^*_j\dfrac { \partial u^*_i(\boldsymbol{x})}{\partial x_j^*} = -\dfrac {\partial p^*}{\partial x^*_i}+\dfrac {\partial }{\partial x_j^*} \bigg ( \frac {1}{\Re } \dfrac {\partial u^*_j(\boldsymbol{x})}{\partial x_j^*} \bigg ),
$$
以及对流-扩散-反应方程（2.7）的控制：
$$
\dfrac { \partial c^*(\boldsymbol{x},t)}{\partial t^*} + \dfrac { \partial c^*(\boldsymbol{x},t) u^*_j(\boldsymbol{x})}{\partial x_j^*} = \dfrac {\partial }{\partial x_j^*} \bigg ( \frac {1}{Pe} \dfrac {\partial c^*(\boldsymbol{x},t)}{\partial x_j^*} \bigg ),
$$
其中$c^*(\boldsymbol{x},t)=c(\boldsymbol{x},t)/c_0$是无量纲浓度，$u^*_j (\boldsymbol{x})= u_{\kern-1pt j} (\boldsymbol{x})/U$是无量纲沿$j$方向的速度，$x_j^*=x_j/d$。我们在容器壁和出口处施加零通量边界条件。雷诺数设定为$\Re = \textit{Ud}/\nu =0.5\lt 1$，并通过确保在$\Re =0.005$时计算的孔隙尺度流动统计（概率分布函数）相同来确认惯性效应可以忽略（此处未显示）。这一结果与Nissan和Berkowitz（2018）的报告一致，他们指出从$\Re =10$开始惯性效应不可忽略。在催化剂表面，用向内指向的法线方向$\lambda$表示，应用了一阶动力学定律（2.8）：
$$
- \frac {\partial c^*(\boldsymbol{x},t)}{\partial \lambda ^*} \bigg |_{S} = {\textit {Da}}_{\textit {I}} \ c^*(\boldsymbol{x},t) \big |_{S},
$$
其中${\textit {Da}}_{\textit {I}}=k_s d/D_m$是孔隙尺度Damk?hler数，它将表面反应的特征速度$k_s$（单位为m s?1）与孔隙尺度扩散速度$D_m/d$联系起来。这里，$\lambda ^*=\lambda /d$表示无量纲界面法线。出于反应器设计的目的，我们定义了一个基于对流的宏观Damk?hler数（2.9）：
$$
{\textit {Da}}_{\textit {II}} = \dfrac {k_v L}{U} = \dfrac {k_s s_p (1-\phi )}{\phi }\dfrac {L}{U},
$$
它将宏观反应速率$k_v = k_s s_p (1-\phi )/\phi$与平均停留时间$L/U$联系起来，其中$s_p=6/d$是催化剂颗粒的比表面积，$\phi \approx 0.5$是孔隙率。我们检查了每种几何配置的两种情况：一个惯性情况，${\textit {Da}}_{\textit {II}} = 0$，以及一个反应性情况，${\textit {Da}}_{\textit {II}} \approx 2$。后者确保反应器在质量控制条件下运行，并实现优化转化（Maggiolo等人2020）。相关无量纲参数总结在表1中，而表2报告了填充分布床反应器应用（如MOST系统）所需的量纲参数。表1. 模拟案例。表格报告了反应器与颗粒大小的比例$R/d$和$L/d$。宏观Damk?hler数设定为${\textit {Da}}_{\textit {II}}=k_v L/U$，分别为0（惯性情况）和2（反应性情况）。所有模拟中的Péclet数为$\textit{Pe} = Ud/D_m \approx 280$，雷诺数为$\textit {Re}=Ud/\nu \approx 0.5$。表2. 表1中模拟的反应器A2和B2的二阶量纲参数，适用于MOST应用等。颗粒直径和反应速率的值反映了典型催化剂的使用情况，容器尺寸和流速与小型微型反应器的尺寸相当（Magson等人2024）。分子扩散率对应于液体中碳氢化合物的典型值（Price和S?derman 2000）。

2.3. 拉格朗日粒子跟踪和沿流线的流体变形
我们应用LPT算法来研究控制流体变形的运动学。实际上，流体变形最终控制了流动传输的标量量的空间分布，即它们的动态混合。从LBM模拟中获得的稳态流场被用作LPT算法的输入，从中提取拉格朗日统计量。LPT的控制方程（2.10）：
$$
\frac {\mathrm{d}\boldsymbol{x}(t,\boldsymbol{X})}{\mathrm{d} t} = \boldsymbol{u}[\boldsymbol{x}(t,\boldsymbol{X}),
$$
使用三阶龙格-库塔方法（Maggiolo、Picano和Guarnieri 2016）求解，其中$\boldsymbol{x}$和$\boldsymbol{X}$分别表示欧拉和拉格朗日位置。我们在入口处随机初始化500条轨迹，并跟踪它们，直到它们平均到达介质长度的一半。提取的拉格朗日统计量用于§ 3.3中量化宏观分散行为，并用于§ 3.2中分析局部流体变形。对于后者，我们采用了Lester等人（Lester、Denz、Le Borgne和de Barros 2018）的方法，该方法涉及构建一个Protean移动参考框架。对于稳态流动，Protean框架本质上是一个流线框架。这种转换使得速度梯度张量变为上三角形状，简化了变形张量的求解，并将拉伸与剪切变形分开。转换表示为两个连续的重定向，$\unicode{x1D64C}(t) = \unicode{x1D64C}_1(t)\unicode{x1D64C}_2(t)$。第一个重定向$\unicode{x1D64C}_1$将速度向量与Protean框架对齐。第二个重定向$\unicode{x1D64C}_2$关于流线方向施加任意旋转。这种旋转角度$\alpha (t)$的特定选择确保速度梯度张量变为上三角形状。角度$\alpha (t)$是从一个常微分方程获得的，其系数取决于速度分量$v_i(t)$（$i=1,2,3$）和在每个时间步骤沿拉格朗日轨迹评估的速度梯度张量分量$\epsilon _{\textit{ij}}(t)$。速度分量可以直接从每个时间步骤的LPT解中获得，而速度梯度张量则通过在计算网格上进行三线性插值来计算。在接近固体表面时，如果插值所需的网格节点位于流体域之外，则通过镜像最近流体节点的速度来施加无滑移边界条件。有关Protean框架的更多细节，请参阅Lester等人的研究（Lester、Denz、Le Borgne和de Barros 2018）。速度梯度张量的推导和验证测试案例在附录A中提供。

3. 结果
3.1. 流动拓扑
我们首先通过表征孔隙尺度的速度和变形来开始分析。图1(b)显示了沿流向$u_x(\boldsymbol{x})$和横向$u_t(\boldsymbol{x})=(u_y^2(\boldsymbol{x})+u_z^2(\boldsymbol{x}))^{1/2}$的孔隙尺度速度的概率密度函数（PDFs）。结果表明，大多数高速流动路径是纵向对齐的。流向和横向速度PDF之间的差异仅弱依赖于所研究的$R/d$比例下的颗粒大小。此外，两种横向约束情况$R/d=4.7$和$R/d=3.1$的速度PDF几乎相同。这表明壁面引起的流体动力效应很小，可以有效地认为介质在横向上是无界的（正如通过壁面到核心的流体比率$d/2/(2R-d)\lt 0.1$（Zió?kowska和Zió?kowski 1988）所确认的）。图2. 通过Delaunay三角剖分计算的孔隙收缩方向$\boldsymbol{\eta }$的极角（$\theta _p$）和方位角（$\varphi _p$）的PDF（图a,b），以及收缩特征尺寸$\sqrt {A_p}$的PDF（图c）。图(d)展示了一个三角形孔隙收缩（顶部），以及相关的孔隙尺度速度$\boldsymbol{u}_p$和法线$\boldsymbol{\eta }$，并显示了用Delaunay三角剖分分割的孔隙空间（底部）。图(a)和(b)中的红色实线对应于$\cos \theta _p$和$\varphi _p$的均匀分布。在(c)中，孔隙收缩的大小显示出狭窄的分布，$\sqrt {A_p}/d \approx 0.5-1$。图(e)展示了流场的横截面，其中无量纲孔隙尺度速度大小$v/U$以对数灰度级别显示：这种表示方法有助于评估孔隙空间的形状，这些孔隙空间在很大程度上是非圆形的。图(f)和(g)分别报告了孔隙空间内椭圆的最小轴与最大轴之比$\varPi_p$的概率密度函数（PDF），以及无量纲最小轴$e_{p,\textit{min}/d$的PDF。这些图表证实了孔隙空间是非圆形的，并且孔隙缩小的概率（$e_{p,\textit{min}$的小值，见图(g)中的红线）与$e_{p,\textit{min}$成线性关系。在最后的图(h)中，我们报告了沿着长度为$L$、直径为$R-r_w$的圆柱体积的平均纵向速度的变化，其中$r_w$是从填充床容器壁的距离。直到$r_w/d \approx 1/2$时，壁效应仍然很明显。在孔隙尺度上，固体障碍物结合无滑移边界条件会产生粘性剪切，这对标量浓度的变形以及调节表面反应非常重要（参考文献：Aquino等人，Ben-Noah, Hidalgo和Dentz2023）。通过应用三维Delaunay三角剖分到粒子中心集上来识别孔隙缩小的特征尺寸。这种方法非常适合球形填充物（参考文献：Ben-Noah, Hidalgo和Dentz，Berkowitz, Cortis, Dentz和Scher2025）。对于由三个粒子中心定义的每个三角形，我们计算了孔隙缩小面积$A_p$、其特征尺寸$\ell_p = \sqrt{A_p}$以及单位法向量$\boldsymbol{\eta}$。后者不一定与正的流动方向对齐。它相对于笛卡尔坐标轴$(x,y,z)$的方向通过极角$\theta_p$和方位角$\varphi_p$来量化，如图2(d)所示。因此，这里使用了三个主要描述符来表征孔隙尺度拓扑的统计特性：（i）孔隙缩小$\ell_p$，（ii）由极角$\theta_p$给出的流动方向对齐，以及（iii）由方位角$\varphi_p$给出的横向方向对齐。图2(a,b)中报告的$\theta_p$和$\varphi_p$的PDF紧密跟随理论上的均匀分布$\cos \theta_p$和$\varphi_p$（红线），表明孔隙法向量没有优选方向。相比之下，图2(c)中的孔隙缩小显示出更狭窄的分布，通常范围从$\ell_p \sim 0.5d$到$d$。这种方向的随机性预计会影响浓度的扩散和混合，这将在§3.6中讨论。

3.2. 拉格朗日流体变形 液体变形由不同的机制引起。虽然剪切流会引起代数变形（参考文献：Souzy等人，Souzy, Zaier, Lhuissier, Le Borgne和Metzger2018），但在随机介质中，拉伸和折叠过程会导致流体元素的指数变形（参考文献：Heyman等人，Heyman, Lester, Turuban, Méheust和Le Borgne2020）。这些机制相互作用，产生变形，表现为幂律函数和指数函数的乘积。为了量化剪切诱导变形与指数变形的相对贡献，我们计算了沿拉格朗日轨迹的变形张量$\unicode{x1D641}^{{\kern1pt} \prime}(t)$。识别主导的变形机制是至关重要的，因为它们控制了溶质的混合和多孔介质的均质化效率。图3. 在体积切片(a)和$(x,z)$平面上的横截面上，欧拉速度大小场$v/U$（b）。(c) 沿拉格朗日轨迹采样的速度$v/U$的PDF。指数$\beta = 3/2$是根据分布的低速尾部使用（3.1）提取的（红线）。深蓝色和浅蓝色圆圈分别对应$R/d=4.7$和$R/d=3.1$的情况，而灰色菱形标记报告了Souzy等人的数据（参考文献：Souzy, Lhuissier, Méheust, Le Borgne和Metzger2020）。我们首先沿流体轨迹提取孔隙尺度速度大小$v=(u_i^2)^{1/2}$，其中$i=1,2,3$表示三个Protean方向。图3(a)和3(b)分别展示了三维切片和对称平面上的欧拉速度大小场。沿拉格朗日轨迹采样无量纲速度$v/U$得到图3(c)中显示的PDF。从这个分布中，我们提取了幂律指数$\beta$，假设低速统计遵循异质介质的特征重尾形式（参考文献：Berkowitz等人，Bijeljic和Blunt2006）\begin{align} P(v\lt U) \propto v^{\beta -1}。\end{align} 这种分析得出$\beta \approx 3/2$（图3 c中的红线）。有趣的是，低速分布并不平坦，如果孔隙截面是圆形的话应该是平坦的（Dentz, Icardi和Hidalgo参考文献：Dentz, Icardi和Hidalgo2018）。在圆形截面上，低速处的PDF遵循管内Poiseuille流的解，其中低速区域（靠近管道壁）比高速区域更可能发生。低速区域的高概率被抛物线流剖面（如Poiseuille流）中低速的低概率所抵消，导致低速的概率分布是平坦的，即$P(v\lt U) \propto v^0$。在非圆形截面上，情况不同：在矩形截面中，随着速度的增加，低速的概率也增加，因为接近墙壁时找到低速区域的概率并不会继续增加（例如参见Stewart和Zhang的参考文献：Stewart和Zhang2013；在二维流动的极限情况下，低速区域的概率是平坦的）。通过在多个截面上进行图像分割（基于$v/U$的阈值分割，见图2 e），我们定量识别和区分了各个孔隙区域；我们观察到孔隙空间主要是非圆形和椭圆形的，其最小轴和最大轴之比$\varPi_p =e_{p,\textit{min}}/e_{p,\textit{max}} \lt 1$，并且孔隙缩小的分布大致遵循$P(e_{p,\textit{min}})\propto e_{p,\textit{min}}$（见图3的(f)和(g)）。de Anna等人（参考文献：Aquino, Le Borgne和Heyman2017）表明，对于二维中孔隙喉部的幂律分布，低速分布遵循的幂律指数的值是孔隙喉部分布指数的一半。我们的发现与这一结果一致，因为（3.1）中的低速幂律指数满足$\beta - 1 = 1/2$。这表明孔隙空间更好地可以用准二维几何形状来近似（其中$e_{p,\textit{min}}$表示二维孔径大小），而不是用具有圆形截面的完全三维孔隙来近似。

图4. (a) 在$R/d=4.7$的情况下，沿拉格朗日轨迹平均的变形张量分量随时间的演变。在早期（$n\lt 10$），纵向剪切诱导的分量$F'_{12}$和$F'_{13}$占主导地位，这一点很好地被预测（3.3）（虚线）所捕捉。插图：计算得到的纵向和横向变形（$\varLambda_\ell$，根据（3.4），$\varLambda_t=\langle \sqrt{{F'_{22}}^2+{F'_{33}}^2+{F'_{23}}^2} \rangle$）与交叉预测（3.6）（点线）的比较。(b) Protean框架（红色箭头）相对于拉格朗日轨迹的方向示例；蓝色阴影区域示意性地表示流体剪切变形$F'_{13}$。(c) 速度梯度张量分量$\epsilon'_{\textit{ij}}(t)$的自相关${c'}^*_\epsilon$，显示了剪切率分量$\epsilon'_{12}$和$\epsilon'_{13}$（纵向剪切）在超过孔隙尺度行程时间$n=1$后的长期相关性，与快速去相关的横向剪切率$\epsilon'_{23}$和指数率$\epsilon'_{22}$、$\epsilon'_{23}$、$\epsilon'_{33}$相比。接下来，我们分析变形张量的分量，其表达式如（3.2）所示。非对角线项$F'_{12}(t)$、$F'_{13}(t)$和$F'_{23}(t)$分别与流动方向（$F'_{12}$、$F'_{13}$）和横向方向（$F'_{23}$）的变形相关，量化了有限时间内的剪切变形以及随后剪切和拉伸的联合效应。对角线项$F'_{11}(t)$、$F'_{22}(t)$和$F'_{33}(t)$捕捉了由速度间歇性引起的沿Protean轴的纯压缩和拉伸。图4(a)显示了它们的时间演变：最初变形主要由剪切控制，$F'_{12}$和$F'_{13}$占主导地位，直到$n \approx 10$（红色阴影区域，$n=tU/d$是特征无量纲时间）。在较晚的时间，指数变形（对角线项$F'_{22}$和$F'_{33}$）在大小上变得相当，并预计将在更长的时间内占主导地位（Lester等人参考文献：Lester, Dentz, Le Borgne和de Barros2018）。在有限系统中，这种剪切机制特别相关，因为宏观长度只有几个粒子直径，可能在指数拉伸发展之前就发生反应。Lester等人的数学框架（参考文献：Lester, Dentz, Le Borgne和de Barros2018）解释了这种行为。在非常早期的时间，变形是弹道性的；在指数拉伸之前（$t\lt t_\epsilon$），剪切占主导。具体来说，$F'_{13}$遵循$n$的幂律（3.3）\begin{align} |F'_{13} (t\lt t_\epsilon )| \approx \gamma ^*_{2} n^{3 - \beta /2}，\end{align}，其中$\gamma ^*_2$是Protean框架$(1,3)$平面中的有效剪切率，$\beta$是速度PDF指数。图4(a)确认了这个模型准确捕获了平均变形$\langle |F'_{13}(t)| \rangle$（其中$\langle \rangle$表示对计算出的拉格朗日轨迹的平均），$\gamma ^*_2 \approx 5$。这个值比有序珠状填充物报告的值（$\gamma ^* \approx 40$）要小（参考文献：Ben-Noah, Hidalgo和Dentz2023）。这种差异是预期中的，因为有序填充物可以触发最佳混合（Turuban等人参考文献：Turuban, Lester, Le Borgne和Méheust2018，Turuban, Lester, Heyman, Borgne和Méheust2019）。$\langle |F'_{12}(t)| \rangle$的演变与$\langle |F'_{13}(t)| \rangle$类似，但略慢，而$|F'_{22}(t)| \ll |F'_{13}(t)|$，这是我们感兴趣的主要结果，表明在有限时间内剪切诱导的变形占主导。代数函数（如剪切变形）在短时间内比指数函数（如拉伸）增长得更快。然而，在长时间内情况则相反。因此，有限时间内从剪切到拉伸的转变是普遍的。哪种变形模式占主导以决定混合取决于具体的剪切强度，我们将在§3.7中进一步讨论。假设$|F'_{12}| \sim |F'_{13}|$，材料线在有限时间内的纵向增长由（3.4）给出\begin{align} \varLambda _\ell (t\lt t_\epsilon ) = \Big \langle \sqrt {{F'_{11}}^2(t)+{F'_{12}}^2(t)+{F'_{13}}^2(t)} \Big \rangle \approx \gamma ^* n^\alpha，\end{align}, 其中$\alpha = 3 - \beta /2 \approx 2.25$，三维中的有效剪切率为$\gamma ^* = (2{\gamma ^*_2}^2)^{1/2} \approx 5\sqrt {2}$。这种缩放将随后用于模拟浓度聚集和混合。

图5. 有限时间Lyapunov指数$\langle \nu (t) \rangle$：通过（3.5）的有限时间解获得的数值估计（标记）与预测（3.7）（$\lambda =0.1$，$\beta =3/2$，红线）的比较（灰色标记：限制更强的情况，$R/d=3.1$）。指数拉伸的开始发生在特征时间$t_\epsilon = 1/\lambda$，其中$\lambda$是无限时间Lyapunov指数（3.5）\begin{align} \lambda = \left \langle \lim _{t \to \infty } \frac {1}{t} \int _0^t \epsilon '_{22}(t),\mathrm{d}t \right \rangle = -\left \langle \lim _{t \to \infty } \frac {1}{t} \int _0^t \epsilon '_{33}(t),\mathrm{d}t \right \rangle。\end{align}。根据Lester等人（参考文献：Lester, Dentz, Le Borgne和de Barros2018），剪切诱导和指数变形之间的转变由（3.6）描述\begin{align} \varLambda _\ell (t) \approx \gamma ^* n^{2 - \beta /2} \exp (\lambda t)，\end{align}，而有限时间Lyapunov指数随时间演变如下（3.7）\begin{align} \langle \nu (t) \rangle \approx \lambda + \bigg (2 - \frac {\beta }{2}\bigg )\frac {\ln (n)}{n}，\end{align}，当$t,n \to \infty$时趋近于$\lambda$。通过（3.5）在有限时间求解$\nu (t)$并用（3.7）进行拟合，得到$\lambda \approx 0.1$（图5，标记和红线）。这个值与随机拉伸方向计算的结果（Lester等人参考文献：Lester, Dentz和Le Borgne2016）相符。来自（3.6）（点线）的预测也与交叉区域观察到的纵向变形（图4插图）相匹配。横向变形$\varLambda _t(t)$应该呈指数级缩放（Lester等人参考文献：Lester, Dentz, Le Borgne和de Barros2018），这一点通过我们的模拟得到了证实（图4 a的插图）。将（3.7）与（3.3）结合，可以直观地观察到有限时间李雅普诺夫指数依赖于主导剪切变形$F'_{13}$（3.8）：
$$
\langle \nu (t) \rangle \approx \lambda + \frac {1}{n} \Big [ \ln \big ( {\gamma ^*_2}^{-1} \big\langle F'_{13}(n) \big\rangle \big ) - \ln (n) \Big ] .
$$
最后，我们注意到中间变形区域主要由相关的剪切事件主导，这些事件在接近指数区域时逐渐失去相关性（见图4c）。特别是，剪切引起的变形在超过孔隙尺度传播时间（$n=1$）后仍然保持相关性。速度梯度分量的自相关性按（3.9）计算：
$$
{c'_{\epsilon }}^* = \frac { \Big\langle \big(\epsilon '_{\textit{ij}}(t)-\langle \epsilon '_{\textit{ij}}(t)\rangle \big)\big(\epsilon '_{\textit{ij}}(0)- \big\langle \epsilon '_{\textit{ij}}(0)\big\rangle \big)\Big\rangle }{\Big\langle \big(\epsilon '_{\textit{ij}}(t)-\langle \epsilon '_{\textit{ij}}(t)\rangle \big)^2 \Big\rangle ^{1/2} \Big\langle \big(\epsilon '_{\textit{ij}}(0)-\langle \epsilon '_{\textit{ij}}(0)\big\rangle \big)^2 \Big\rangle ^{1/2},
$$
其中$\epsilon '_{\textit{ij}}(t)$是上三角速度梯度张量的分量。这种相关性在控制标量聚集的随机性中起着关键作用（参见Duplat等人，参考文献Duplat, Jouary和Villermaux2010），并将在§3.6中进一步讨论。

3.3. 示踪剂的扩散和浓度前沿的演变
我们现在转向标量扩散和混合的分析。我们首先分析惰性多孔介质，然后再将讨论扩展到反应性系统。在本节和§3.4中，我们根据孔隙尺度模拟的欧拉数据来分析扩散和混合，而在§3.5、3.6和3.7中，我们使用基于拉格朗日变形的混合模型来解释结果。图6显示了惰性多孔反应器中标量传输的快照（$R/d=4.7$）。无量纲浓度$c^*(\boldsymbol{x},t)=c(\boldsymbol{x},t)/c_0$以恒定浓度$c_0$在左边界注入。在平均前沿位置$c^*\approx 0.5$（蓝色区域）附近可以看到浓度丝状物（或片状物）。当一个标量被注入多孔反应器时，流动变形强烈影响浓度的分布。这种效应在平均浓度前沿附近尤为明显，即靠近等浓度面$c^*\approx 0.5$（见图6）。这里，浓度丝状物沿着高速度流动路径形成，这些路径对应于变形强烈的区域。更具体地说，浓度片状物优先沿着主导流动方向形成和拉伸。对于补充浓度场$c_w=1-c^*$也会发生同样的现象，我们将其称为白色浓度元素（与黑色高浓度区域相对）。最显著的白色区域出现在催化剂粒子的尾流中（见图7a）。在横截面上，这些区域类似于二维片状物（见图7b），它们偶尔会横向重叠（红色箭头）。

图7. (a) 显示了以平均前沿位置$Ut/d=n$为中心的扩散体积的纵向延伸$n^\delta$的快照。低浓度区域（白色）形成在粒子尾流中（红色箭头）。(b) 横截面显示了白色区域的片状三维结构（面板a中的红色箭头）。(c) 放大了面板(b)中红色矩形内的两个重叠片状物。这里的白色片状物宽度为$s_w$，显示为灰色。(d) 观察到的片状物是由几个基本浓度片状物组成的束，其最大值为$c_e$，厚度为$s_e$。这些束的平均厚度为$s_w$，浓度为$c_w$。随着$n = Ut/d$，平均浓度前沿向前推进，这也代表了遇到的固体障碍物的平均数量。白色浓度元素在平均前沿附近积聚，随着时间的推移而延长，并从后方增加到头部。这导致初始步骤剖面（$c^*=1$，$c_w=0$在后方；$c^*=0$，$c_w=1$在头部）逐渐转变为更平滑的纵向剖面。这种平滑是由于纵向扩散造成的，它将浓度元素沿着流向方向扩散。为了量化浓度前沿的演变，我们从拉格朗日示踪剂的平均平方位移（MSD）计算纵向扩散（3.10）：
$$
\textit {MSD}_\ell = \dfrac {\big\langle [ (x_k(t)-x_k(0))-\langle x_k(t)-x_k(0)\rangle ]^2 \big\rangle }{d^2},
$$
（3.11）
$$
\textit {MSD}_t = \dfrac {1}{2d^2} \Big ( {\big\langle [ (y_k(t)-y_k(0))-\langle y_k(t)-y_k(0)\rangle ]^2 \big\rangle } + {\big\langle [ (z_k(t)-z_k(0))-\langle z_k(t)-z_k(0)\rangle ]^2 \big ),
$$
其中$\textit {MSD}_\ell$和$\textit {MSD}_t$分别表示纵向和横向的MSD，$(x_k(t),y_k(t),z_k(t)$是拉格朗日示踪剂的位置。特征流动时间$n=1$后，纵向MSD表现出超扩散尺度，而横向MSD是扩散的（见图8）。

图8. 从拉格朗日轨迹计算出的平均平方位移（3.11）。对于(a) $R/d=4.7$和(b) $R/d=3.1$，纵向MSD表现出超扩散行为。(c) 比较了两种计算的扩散长度：(i) 从拉格朗日统计得出的$\sqrt {\textit {MSD}_\ell }$（蓝色虚线），以及(ii) 从欧拉场得出的$2(n-x_1)$（标记），其中$n$和$x_1$是平均值和最低前沿位置（见图7a）。纵向扩散将示踪剂沿着平均流动方向扩散，产生的扩散体积的流向长度也可以从欧拉浓度场估计出来。我们将这个体积的最低位置$x_1$确定为$c_w$的平均值严格保持正值的第一个横截面。然后，我们将从欧拉浓度场估计的纵向扩散长度$2(n-x_1)$（见图7a中的示意图）与拉格朗日评估$\sqrt {\textit {MSD}_\ell } \propto n^\delta$进行比较，其中$\delta ={3/4}$是从图8中显示的超扩散趋势中提取的。这两种测量方法相当吻合（见图8c），表明模拟的欧拉数据和拉格朗日数据之间的一致性。因此，纵向扩散系数与平均速度成比例，$D_\ell \approx Ud$，浓度前沿在时间上呈准线性扩展的体积内扩散，其长度随$n^\delta$增加，其中$\delta \approx 0.75$。在更长的时间后，我们预期这种超扩散状态将转变为菲克扩散，这与各向异性（Maggiolo等人，参考文献Maggiolo, Picano和Guarnieri2016）和高度异质性介质（Souzy等人，参考文献Souzy, Lhuissier, Méheust, Le Borgne和Metzger2020）中的观察结果一致。线性$D_\ell - U$比例关系在$\textit{Pe}=O(10^2)$时是合理的，这与之前的研究结果一致，其中纵向扩散系数预计会按$D_\ell \propto Pe^{1.19}$比例变化，直到$\textit{Pe} \approx 400$，然后在更高的$\textit{Pe}$时转变为线性比例（Bijeljic, Muggeridge和Blunt，参考文献Chiogna, Hochstetler, Bellin, Kitanidis和Rolle2004）。横向扩散应显示出更稳健的关系，按$D_t \propto Pe$比例变化（Bijeljic和Blunt，参考文献Bijeljic, Muggeridge和Blunt2007）；在随机堆积床中，这种行为已经实验证明在$\textit{Pe} \gtrapprox 10^3$时观察到（Scheven，参考文献Scheven2013）。值得注意的是，作为另一种方法，也可以通过连续时间随机游走理论从速度分布计算示踪剂的纵向和横向扩散，该理论提供了沿流线的欧拉和拉格朗日速度之间的关系。更多阅读材料，请参阅Dentz等人（参考文献Dentz, Kang, Comolli, Le Borgne和Lester2016）。

3.4. 孔隙尺度混合的演变
为了表征混合，我们使用浓度方差作为混合动力学的度量。在孔隙尺度上，其他常见的度量包括稀释指数（Chiogna等人，参考文献Dentz和de Barros2012）。如果观测分辨率足够高，浓度方差能够捕捉到小尺度变化。这可以在固定的或移动的纵向区域（“流体支撑”）计算，或者在该区域内浓度超过阈值的子集内计算（“浓度支撑”）（Lester等人，参考文献Lester, Dentz和Le Borgne2016）。在这项工作中，我们采用流体支撑的方差作为主要的混合度量。我们关注一个靠近平均浓度前沿的子体积，定义为中心浓度$c^*=0.5$，该子体积沿着反应器移动。实际上，我们选择一个厚度为$md$的离散圆盘（其中$m$是一个常数，介于4到6之间，$d$是粒子大小），使得体积平均浓度为$\mu _v = 0.5$。这种选择减少了计算方差时平均时间波动的影响。子体积大约以$n = Ut/d$的速度移动。为了避免明显的壁效应（也见图2的面板（h）），我们将圆盘的面积$A$限制为半径为$R-d/2$的圆。为了清晰起见，我们定义了一个移位的流向位置$x'' = (n-m/2)d$（见图9右面板），这是位于平均前沿位置$n$的子体积内的长度为$m$的位置。然后计算子体积内的体积方差（3.12）：
$$
\sigma ^2_v(t) = \dfrac {1}{Amd} \int _{0}^{md} \int _A {c}^2(\boldsymbol{x}_f,t) \mathrm{d}A \ \mathrm{d}x'' \ - \mu _v^2(t) ,
$$
其中$c(\boldsymbol{x}_f,t)$是流体位置$\boldsymbol{x}_f$的浓度，$\mu _v$是其体积平均浓度。为了分离纵向扩散和孔隙尺度混合，我们引入空间平均运算符$\langle \ \rangle _{\kern-1pt j}$，应用于厚度为$d$、面积为$A$、中心位于$x''/d = (2j-1)/2$、$j=1,\ldots ,m$的三维圆盘。然后（3.12）中的第一和第二矩变为（3.13）：
$$
\mu _v = \dfrac {1}{m} \sum _{j=1}^m \dfrac {1}{Ad} \int _{(j-1)d}^{jd} \int _A c(\boldsymbol{x}_f,t) \mathrm{d}A \ \mathrm{d}x'' = \dfrac {1}{m} \sum _{j=1}^m \langle c \rangle _{\kern-1pt j}
$$
$$
\dfrac {1}{Amd} \int _0^{md} \int _A {c}^2(\boldsymbol{x}_f,t) \mathrm{d}A \ \mathrm{d}x'' = \dfrac {1}{m} \sum _{j=1}^m \big\langle c^2 \big\rangle _{\kern-1pt j} .
$$
将（3.12）和（3.13）结合，浓度方差可以重写为（3.14）：
$$
\sigma ^2_v(t) = \dfrac {1}{m^2} \sum _{j=1}^m \sum _{i=1}^{j-1} \big ( \langle c \rangle _{\kern-1pt j} - \langle c \rangle _i \big )^2 + \dfrac {1}{m} \sum _{j=1}^m \big\langle {c'}^2 \big\rangle _{\kern-1pt j} ,
$$
其中$c' = c - \langle c \rangle _{\kern-1pt j}$。右侧的第一项量化了一维纵向扩散，测量了横截面平均浓度之间的差异。它与纵向浓度梯度的平方$(c_0^*/L_\ell )^2$成比例，其中$L_\ell ^2 \sim \textit {MSD}_\ell \propto n^{2\delta }$（§ 3.3），因此对于$\delta \sim 0.75$，其比例为$n^{-2\delta } \approx n^{-1.5}$。第二项捕捉了横向变化，测量了孔隙尺度差异和混合，通过计算$\langle \boldsymbol{\cdot }\rangle _{\kern-1pt j}$在厚度为$\lessapprox d$的圆盘上消除了纵向孔隙间的差异，$\lessapprox d$是速度波动的特征长度尺度。

图9. 对浓度方差$\sigma _v^2$的贡献（3.14）对于(a) $R/d=4.7$和(b) $R/d=3.1$（平均两个实现）。纵向扩散（空心方块）随时间以$n^{-2\delta }$减少，而孔隙尺度混合（实心圆）衰减得更慢$\sim n^{-1/2}$（实线蓝色，见（3.27））并在长时间内占主导地位（内嵌图）。面板(c)展示了用于计算（3.14）的三维圆盘的划分。图9显示了两种贡献对方差的影响。纵向扩散减少了沿流动方向的浓度差异，按$n^{-2\delta }$比例缩放（空心方块）。孔隙尺度混合衰减得更慢（实心圆）并在子体积中占主导地位。这种行为在有限时间内很重要，将在后续章节中进一步分析。

3.5. 剪切流中的层状模型
我们观察到浓度区域沿流向方向明显延长，现在我们探讨它们是如何在空间中形成和排列的。为此，我们检查了白色浓度的分布。这些低浓度区域可以被视为几个基本浓度片状的束或聚集体，这些片状物起源于粒子尾流。我们用$c_e$表示基本片状的最大浓度，用$c_w$表示束内平均厚度为$s_w$的最大浓度。图7的面板(c)和(d)展示了一个束形成的草图。每个基本片状物起源于粒子尾流，并随后被流动变形。我们可以使用层状模型（Villermaux，参考文献Villermaux2019；Dentz等人，参考文献Dentz, Hidalgo和Lester2023）来模拟它们在有限时间内的演化，即在主导变形机制下，即剪切流中的演化。剪切引起的延长按$\varLambda _\ell \propto n^\alpha$比例变化，其中$\alpha$是从（3.4）中的拉格朗日运动学中提取的。每个片状物在短混合时间$t_m$后的最大浓度减少为$c_e = c_0/(s_e n^\alpha )$，其中$s_e\propto n^{1/2}$是横向尺寸。由于 $\alpha \approx 2.25$ 大于色散体积增长的指数 $\delta \approx 0.75$（$n^\alpha \gt n^\delta$），片层会因为剪切变形而重叠，其速率为 $n^{\alpha -\delta +1/2}$，从而产生有效浓度 $c_e n^{\alpha -\delta } s_e$。由于它们的密度逐渐增加，重叠现象更加明显：由于它们沿着流动方向延伸，每个截面的基本片层数量随之增加，即 $n_e = x' n_p$，从稀疏区域（$x'\sim 0$）变为密集区域（$x'\sim n$）。我们强调浓度片层的密度会随时间增长，因为更多的片层到达给定的纵向位置。基本片层的聚集从而形成了束状结构（见图10c和10f），并且每个 $x'$ 处的质量守恒需要满足方程（3.15）：
\begin{align}
\langle c_{w} \rangle \langle s_w \rangle \dfrac {n_w}{n_p} = c_e n^{\alpha -\delta } s_e \dfrac {n_e}{n_p} = \dfrac {x'}{n^\delta } ,
\end{align}
其中 $c_0=1$，符号 $\langle \rangle$ 表示每个截面的样本平均值。我们通过孔隙尺度浓度场对这一模型进行数值验证。我们追踪位于区间 $x'_1=[x_1 : x_1+ (n d))^\delta ]$ 内不同截面的最大值 $c_w$，该区间对应于以平均前沿位置 $Ut/d=n$ 为中心的纵向色散长度（见图7a）。图10的上半部分显示了不同无量纲时间 $n$ 下最大值的样本平均值 $\langle c_w \rangle$，这个值随着纵向位置 $x' = x'_1-x_1$ 在色散体积内的增加而增加。图10：(a,b) 上半部分：不同时间 $n$ 下沿色散体积纵向发展的白色束状结构中最大浓度的样本平均值 $\langle c_{w} \rangle$，对于 $R=4.7$（a）和 $R=3.1$（b）；中间部分：平均横向束状结构尺寸 $s_w$，它随着 $n$ 的增加而增加（蓝线）；下半部分：质量守恒（3.15）的验证（红线）；(c) 低浓度 $c_e$ 的基本片层示意图，这些片层被流动纵向拉伸，并重叠形成浓度为 $c_w$、厚度为 $s_w$ 的束状结构；(d) 束状结构嵌入在长度约为 $n^\delta$ 的色散体积中（蓝色阴影区域）。在平均前沿位置的上下游，白色浓度的稀疏和密集区域被标出（红色垂直线）；(e) 当束状结构的宽度 $s_w$ 超过它们之间的间距 $\ell _w$ 时，会发生重叠，产生一个新的最大值（红点），此时 $n_w/n_p \sim s_w/\ell _w$。在最右侧的面板（f）中，两个浓度场的快照显示了白色浓度基本片层在下游位置的出现，先是形成束状结构，随后数量增加并合并（底部）。在稀疏区域，束状结构由少量片层组成，聚集程度较低，导致 $\langle c_w \rangle$ 随时间的延长和扩散而呈幂律减小（见图10上半部分）。在密集区域，聚集占主导地位，重叠使得 $\langle c_w \rangle \to 1$。扩散作用横向稀释了白色束状结构。平均横向扩散距离 $\langle s_w \rangle$ 可以通过每个孔隙的最大值数量 $n_w/n_p$ 和它们的平均间距 $\langle \ell _w \rangle$ 来推断，其中 $n_p=(2R/d)^2$ 是每个截面的平均孔隙数。局部最大值出现在不同束状结构的交叉点，预计随着束状结构之间的重叠增加，其数量会随之增加，估计为 $\sim \langle s_w \rangle /\langle \ell _w \rangle$。这导致了方程（3.16）：
\begin{align}
\langle s_w \rangle \sim \langle \ell _w \rangle \dfrac {n_w}{n_p} \propto n^{1/2} ,
\end{align}
，这意味着如果基本片层呈扩散性增长（$s_e \propto t^{1/2}$），则浓度会随时间扩散性增长。（在（3.16）中，$\ell _w$ 和 $s_w$ 都被粒径 $d$ 标准化了。）从图10的中间部分我们可以观察到，在固定 $x'$ 的条件下，$\langle s_w \rangle$ 随 $n$ 的增加而增加，这支持了基本片层扩散性增长以及随后束状结构重叠的解释。最后，在图10(a)和10(b)的最下方面板中，我们也验证了方程（3.15）。鉴于模拟结果与剪切流中层状变形模型的良好吻合，接下来我们将更详细地分析混合动力学，并通过层状聚集提供更精确的预测。特别是，我们关注浓度最大值的第二矩，它提供了孔隙尺度混合的准确度量。

3.6. 剪切诱导变形下的聚集模型
我们采用剪切流变形模型（Souzy等人，参考文献Souzy, Zaier, Lhuissier, Le Borgne和Metzger2018；Villermaux，参考文献Villermaux2019）来描述有限时间内的混合动力学。我们之前指出，该模型与白色浓度片层的第一矩是一致的：正是剪切流和扩散共同作用，在粒子尾流处诱导了白色浓度片层的形成，并在有限时间内实现了连续的混合。这里我们简要回顾一下剪切流变形模型。在非常早期的时间，平流作用占主导，流动引起的拉伸使得浓度片层沿横向方向压缩。由于剪切诱导的变形占主导（见§3.2），这些白色片层的横向厚度 $s_e(t)$ 根据方程（3.17）减小：
\begin{align}
s^{(1)}_e (t) \approx s_0 (\gamma t)^{-\alpha } ,
\end{align}
其中 $s_0$ 和 $s_0^*=s_0/d$ 分别表示有量纲和无量纲的初始厚度。图11：(a) 在初始时刻，由于剪切流，低浓度片层从粒子尾流中形成。它们的初始横向尺寸为 $s_0=d/2$，因此总浓度 $c^*=0.5$；(b) 孔隙尺度混合的示意图：在混合时间 $t_m$，片层达到了它们的最小厚度 $s_b$；混合时间 $t_m$ 之后，扩散作用使它们合并，形成最大浓度为 $c_w$、厚度为 $s_w$ 的束状结构（c）。在(d)中，两种重叠模式被示意出来：浓度片层I由于顺流方向的拉伸而与片层II纵向重叠，而片层II和III之间的横向重叠则由扩散维持。应用Ranz的变换（Ranz参考文献Ranz1979）并在变换后的拉格朗日框架中求解相关的平流-扩散方程，可以估计特征混合时间 $t_m$，在这个时间点扩散（$D/s^2(t)$）开始抵消平流引起的压缩（$-(1/s_e)(\mathrm{d}s_e/\mathrm{d}t)$（3.18）：
\begin{align}
t_m \sim \dfrac {1}{\gamma } \big(\gamma ^* {s^*_0}^2\textit {Pe} \big)^{1/(2\alpha +1)} ,
\end{align}
其中 $\gamma = \gamma ^* U/d$ 是三维空间中的剪切诱导变形量，如§3.2所定义。使用测量的 $\gamma ^*= 5 \sqrt {2}$ 和 $\alpha =2.25$，我们得到 $t_m \approx d/(3U)$。在这个混合时间 $t_m$，拉伸和压缩作用产生了具有横向厚度 $s^{(1)}_b \approx s_0 \big(\gamma ^* {s_0^*}^2 \textit {Pe}\big)^{-\alpha /(2\alpha +1)}$ 的白色浓度片层，随后这些片层通过扩散作用扩展为 $s_e^{(1)} \sim \sqrt {2Dt}$。这些浓度片层主要沿流线方向延伸，因此与流线方向对齐（见图11，面板a）。在 $t_m$ 之前，靠近前沿的区域，拉伸作用将来自不同流向位置的越来越多片层聚集在一起，而压缩作用阻止了它们的重叠。超过 $t_m$ 之后，横向扩散允许这些元素以 $t^{1/2}$ 的速率扩展，它们的最大浓度根据方程（3.20）衰减：
\begin{align}
c_e \approx \dfrac {c_0 \ s_0}{(\gamma t)^\alpha \sqrt {2Dt}}\propto n^{-(\alpha +1/2)} .
\end{align}
在这个阶段，由于剪切诱导的变形和横向扩散的共同作用，这些元素可以自由重叠，如图11所示。早期变形的关键积极效果在于它增强了混合：通过横向压缩浓度元素，平流作用使它们更加接近，从而更有效地促进了扩散作用下的合并。因此，我们基于剪切诱导的拉伸作用定义了一个早期聚集模型，这是有限时间内的主导变形机制，而在更长时间内，指数拉伸占主导。在剪切诱导的变形下，当 $t \gt t_m$ 时，白色浓度片层的聚集变得显著，这是由于两个主要机制：（i）片层的纵向生长与色散体积之间的不匹配，发生率为 $n^{\alpha -\delta }$；（ii）片层密度的增加和它们物理间距的减小，发生率为 $N(t)$。物理上的接近性增加导致了重叠，因为更多的纵向延伸的浓度元素占据了相同的孔隙空间（见图6），而横向扩散作用扩大了它们的有效厚度。每个孔隙中的重叠总数根据方程（3.21）估计为：
\begin{align}
N(t\gt t_m)\approx (\gamma t)^\alpha \sqrt {2Dt}/d \propto n^{\alpha +1/2} ,
\end{align}
这主导了聚集动力学，因为 $N \gt n^{\alpha -\delta }$。因此，聚集体的平均浓度 $c_e N$ 在前沿附近保持大致恒定。质量守恒也反映在白色束状结构最大值的第一矩 $\langle c_w(x=n) \rangle$ 中，如图10所示。尽管如此，这种恒定性不适用于方差。鉴于前沿附近的平均浓度是恒定的，可以应用Duplat和Villermaux的随机聚集模型（参考文献Duplat和Villermaux2008）来估计方差的衰减，该模型基于浓度片层之间重叠事件的计数。根据这个模型，可以通过随机添加独立且分布相同的变量，或者在它们相关时通过均值和方差的总和来模拟不相关浓度片层之间的聚集（Duplat等人，参考文献Duplat, Jouary和Villermaux2010）。在前一种情况下，利用中心极限定理，基于假设孤立片层或束状结构内的浓度是均匀分布的，并在执行聚集时随机选择（也见Weiss和Provenzale，参考文献Weiss和Provenzale2007）。在 $t_m$ 时每个孔隙中的片层数量是它们早期延伸超过孔隙空间的结果（3.22）：
\begin{align}
n_b = (\gamma t_m)^\alpha ,
\end{align}
结合（3.18）和（3.19），定义了片层之间的平均无量纲间距（3.23）：
\begin{align}
\ell _b = 1/n_b = s_b/s_0 .
\end{align}
在 $t_m$ 之后，片层通过两种机制发生重叠：持续的纵向剪切诱导拉伸（$N_\ell$）和横向扩散（$N_t$），因此总的重叠数为 $N = N_\ell N_t$（见（3.21）（3.24）：
\begin{eqnarray}
&&N_\ell (t\gt t_m) = (\gamma t)^\alpha \ell _b ,
\nonumber \\
[5pt]
&&N_t(t\gt t_m) = s_e/(\ell _b d) \approx \sqrt {2Dt}/({\rm d}\ell _b) .
\end{eqnarray}
单个孤立浓度元素的方差与其最大浓度的平方成正比（Meunier和Villermaux，参考文献Meunier和Villermaux2010）（3.25）：
\begin{align}
\sigma ^2_e(t) \approx c_e^2 \propto n^{-2\alpha -1} .
\end{align}
纵向重叠将这些元素聚集成束状结构（见图10，面板c和f）。由于剪切诱导的变形在大于孔隙尺度 $d$ 的长度上是相关的（§3.2），束状结构的方差可以近似为 $N_\ell$ 个单独方差的相关和（3.26）：
\begin{align}
\sigma _w^2(t) \approx c_{w,\textit {max}}^2 = N_\ell ^2 \sigma _e^2 \propto n^{-1} .
\end{align}
图12：(a) 模拟结果（黑色标记：$R/d=4.7$；蓝色标记：$R/d=3.1$）与随机横向聚集的预测（（3.27）以及指数拉伸模型（虚线）的比较。另一个具有双反应器长度的模拟结果也显示出来（空心方块标记），然而直到 $n \approx 50$ 时几乎观察不到指数衰减。（b）根据方程（3.17）（蓝色实线）和指数拉伸（红色虚线）预测的浓度片层压缩，标记显示了混合时间 $s_b^{(1)}$（蓝色）和 $s_b^{(2)}$（红色）时达到的厚度，对于 $\textit {Pe} = 3 \times 10^2$ 和 $3 \times 10^4$ 分别。（c）验证了体积方差和孔隙尺度混合项 $\sim \langle c_{w,\textit {max}}^2 \rangle$。相比之下，横向重叠在截面上是随机发生的。横向重叠产生的浓度值是属于不同束状结构的浓度之和，这些束状结构在截面空间中有随机的方向。这种随机聚集规则可以通过§3.1中观察到的方位孔隙口的随机方向来解释，并在图2中报告。根据这种推理，前沿附近的最大浓度是由根据束状结构的特征PDF分布的 $N_t$ 个浓度值的随机叠加产生的。利用(3.20)和(3.24)，我们得到(3.27)：
\begin{align}
\sigma ^2(t) \approx N_t \sigma _w^2 = N_t N_\ell ^2 c_e^2 = B n^{-1/2},
\end{align}
其中$N_t(t)$表示扩散重叠数量的演变，$B$是一个常数，取决于初始条件$c_0, s_0$、剪切引起的变形率$\gamma ^*$和Péclet数，如(3.28)所示：
\begin{align}
B = \dfrac {c_0^2 s_0}{\sqrt {2}} \Bigg ( \dfrac { {s_0^*}^2Pe}{{\gamma ^*}^{2\alpha }} \Bigg ) ^{1/(4\alpha +2)} .
\end{align}
为了更好地突出Pe数的依赖性，可以将(3.27)重写为(3.29)：
\begin{align}
\sigma ^2(t) \propto ({\gamma ^* Pe})^{-\alpha /(2\alpha +1)} t^{-1/2},
\end{align}
在线性剪切的情况下（$\alpha =1$），这产生了与混沌重启模型（Aquino等人，参考文献Ben-Noah, Hidalgo和Dentz2023）预测的物种到达多孔介质固体表面的旅行时间方差相同的缩放关系。当前的分析表明，可以扩展这些模型以考虑非线性剪切和$\alpha \ne 1$的情况。数值结果证实，前沿附近的方差可以很好地用束状最大值的方差来近似，即$\sigma ^2(t) \sim \langle c_w^2\rangle$（见图12，c面板）。使用从拉格朗日变形中提取的$\gamma ^*$值（§ 3.2）计算(3.28)得到$B \approx 0.2$，将其代入(3.27)后，得到的预测与孔隙尺度混合的欧拉数据（图9，蓝色实线对比实心标记，大面板和插图）非常吻合。这证实了有限时间内的混合主要由剪切引起的变形所主导。所提出的聚类模型结合了相关的纵向和随机的横向重叠，预测浓度方差的衰减比完全三维随机聚类要慢，后者的方差按$n^{-3/2}$缩放。然而，早期的拉伸和压缩加速了横向重叠，与纯扩散相比，导致方差衰减与平均速度而不是分子扩散系数$D$成正比，这与$n \propto U$一致。

3.7. 剪切引起的变形与指数拉伸
在更长的时间范围内，指数拉伸开始占主导地位。然而，在指数拉伸的特征开始时间之后，即$n \sim 1/\lambda \approx 10$，缓慢的剪切引起的混合仍然主导着整个混合过程。在指数拉伸期间，浓度元素必须达到一个与初始浓度无关的稳定最小厚度，Villermaux（参考文献Villermaux2019）：
\begin{align}
s^{(2)}_b \approx d (\lambda \textit {Pe})^{-1/2}.
\end{align}
将无量纲剪切引起的变形率$\gamma ^* = \gamma d / U \approx 14$、$\alpha \approx 2.25$和Lyapunov指数$\lambda = 0.1$的值代入(3.19)和(3.30)，我们发现$s^{(1)}_b \lt s^{(2)}_b$，因为(3.31)：
\begin{align}
\lambda \ll {\gamma ^*}(\gamma ^*{s_0^*}^2 \textit {Pe})^{-1/(2\alpha + 1)}.
\end{align}
特别是，我们发现$s^{(2)}_b \approx 5 s^{(1)}_b$。这种尺度分离$s^{(1)}_b \lt s^{(2)}_b$适用于中等的Péclet数，表明浓度元素首先通过剪切引起的变形达到$s_b^{(1)}$，然后在后来的指数拉伸下通过扩散聚集，最终达到$s_b^{(2)}$的尺度（前提是剪切引起的混合时间$t_m$，这里是$t_m=0.33 d/U$，小于指数拉伸的开始时间$t_\epsilon =d/(U\lambda )$）。在受到指数拉伸之前，浓度场通过扩散重新组织。为了研究预期的渐近行为，我们通过将反应器长度加倍到$L/d = 84$进行了额外的模拟，这允许计算方差衰减直到$n \approx 50$。结果显示在图12中。观察到方差衰减有轻微的变化，暗示着向指数拉伸的过渡。然而，有限的计算域阻止了明确的确认。对于紧凑的、有限大小的多孔反应器，其中反应通常发生在$n \approx 10$ – $50$范围内，我们得出剪切引起的变形仍然是主导的混合机制。将早期观察到的方差的代数衰减与类似随机堆积床中的实验观察结果进行比较是有启发性的（Heyman等人，参考文献Heyman, Lester, Turuban, Méheust和Le Borgne2020，参考文献Heyman, Lester和Le Borgne2021），其中混沌流动产生了方差的指数衰减。Lyapunov指数（$\lambda \approx 0.17$）和混沌混合的开始时间（$t_m^{(2)} \sim 6$）与我们的发现质量上一致。然而，在这些实验中没有观察到早期的代数衰减。这种差异可以归因于不同的注入协议。我们还注意到，由于那些实验中的Péclet数比我们的模拟高出一到两个数量级，因此(3.31)预测的尺度分离大大减小，使得任何可能的有限时间剪切变形与指数变形相比不那么明显。为了更好地突出$\textit{Pe}$的依赖性，我们可以明确地将(3.31)用$\textit{Pe}$表示。然后，为了在$s^{(1)}_b$和$s^{(2)}_b$之间获得相关差异，例如在当前模拟中$s^{(2)}_b \gt 4 s^{(1)}_b$，Péclet数必须满足条件(3.32)：
\begin{align}
\textit {Pe} \lt \frac {{\gamma ^*}^{2\alpha }}{ {s_0^*}^2 (4^2 \lambda )^{2\alpha +1}},
\end{align}
这个条件取决于剪切指数$\alpha$。例如，以特征孔径$s_0^*=1$作为初始浓度厚度，$\lambda =0.1$和$\gamma ^*\approx 10$，对于强非线性剪切（$\alpha =2$），我们发现$\textit{Pe} \lt 10^3$，而对于线性剪切（$\alpha =1$），我们有$\textit{Pe} \lt 2 \ 10^1$。方程(3.32)可以解释如下：如果Péclet数远小于表征指数混合和剪切引起的变形开始时间尺度的比率，剪切在有限时间内主导混合。

3.8. 带有表面反应的混合
在存在反应性表面的情况下，例如流体-催化剂颗粒界面（$\textit {Da}_{\textit {II}} \gt 0$），表面浓度的均匀化和局部反应由横向重叠驱动的混合机制控制，该机制决定了(3.27)中体积方差的缩放。实际上，我们观察到，对于惰性和反应性介质，流体-固体界面的浓度方差$\sigma _s^2$，通过平均值$\mu _s$归一化，遵循(3.33)：
\begin{align}
\sigma _s^2/\mu _s^2 \ (t\lt t_k) = A n^{-1/2}/\mu _s^2(0) \propto ({\gamma ^* Pe})^{-\alpha /(2\alpha +1)} (t_D/t)^{1/2},
\end{align}
直到特征反应时间$t_k = 1/k_v$（其中$t_D=d^2/D$是一个特征扩散时间）。在反应性情况下，我们测量到略微更高的值，但时间趋势保持不变。在更晚的时间（$t \gt t_k$），对于反应性情况，$\sigma _s^2 / \mu _s^2$似乎饱和（见图13下方的实心标记）。这种效应是由反应性表面处的局部浓度耗尽所决定的，这抵消了混合作用。然后，假设这种机制在长时间内抑制了反应性表面的混合，我们可以在时间$t_k$使用(3.33)来估计长时间内的方差-平均值比率(3.34)：
\begin{align}
\sigma ^2_s/\mu ^2_s (t\gt t_k) \propto ({\gamma ^* Pe})^{-\alpha /(2\alpha +1)} {{\textit {Da}}_{\textit {I}}}^{1/2},
\end{align}
这再次与混沌重启模型中的Pe缩放一致（Aquino等人，参考文献Ben-Noah, Hidalgo和Dentz2023），前提是设置$\alpha =1$。这种一致性表明，随机聚类模型和混沌重启模型原则上可以捕捉到反应性多孔介质中流体-固体表面局部反应性的异质分布。图13的(c)面板展示了对于粒子大小$R/d = 3.1$和$R/d = 4.7$的比较，分别对应于$L/d = 28$和$L/d = 42$（见表1），其中$\textit {Da}_{\textit {I}} = 3.33$和$2.22$。我们注意到较低的$\textit {Da}_{\textit {I}}$数（较慢的反应）允许更多的混合发生。出于设计目的，值得注意的是$\sigma ^2_s \propto \mu ^2_s {{\textit {Da}}_{\textit {II}}}^{1/2} (d/L)^{1/2}$，并且表面反应的均匀化（低$\sigma _s^2$）随着系统长度与粒子大小的比率$L/d$的增加而增加。因此，对于给定的催化剂反应速率$k_v$，最佳的反应器设计应该（i）通过保持$\textit {Da}_{\textit {II}} = k_v L / U = O(1)$来实现接近完全转化，以及（ii）通过减小粒子大小$d$或按比例增加流速$U$和反应器长度$L$来促进有效混合（低$\sigma _s$）。实际上，由于生产限制和孔堵塞的风险，减小粒子大小往往是具有挑战性的，这限制了反应器的效率。另一方面，增加流速和反应器长度提供了一种更实际的策略——我们最近在微反应器实验中成功验证了这种策略（Magson等人，参考文献Magson, Maggiolo, Kalagasidis, Henninger, Munz, Kn?bbeler-Bu?, H?lzel, Moth-Poulsen, Funes-Ardoiz和Sampedro2024）。图13显示了（a）$R/d=4.7$和（b）$R/d=3.1$的惰性（顶部面板）和反应性（下部面板）介质中的流体-固体界面浓度平均值$\mu _s$和方差$\sigma _s^2$的演变。比率$\sigma _s^2/\mu _s^2$的衰减相似，遵循(3.27)（实线蓝色曲线），在惰性和反应性情况下直到特征反应时间$t_k=1/k$（底部面板中的垂直红色虚线）。在$t_k$之后，反应性位置的方差-平均值比率$\sigma _s^2/\mu _s^2$饱和。在(c)面板中比较了两个$L/d$情况下的$\sigma _s^2/\mu _s^2$，表明在相同的$\textit {Da}_{\textit {II}}$条件下（也见表1），较长反应器的混合更有效，正如(3.34)所预测的。

4. 结论
在这项研究中，我们通过结合基于格子-玻尔兹曼求解器的拉格朗日和欧拉分析，研究了有限时间内多孔介质中的标量传输和混合，揭示了流动引起的变形、扩散和反应过程之间的相互作用。尽管混合是混沌的，但剪切引起的变形可以通过横向压缩浓度区域和促进纵向聚集来主导早期混合动态，从而增强混合。横向扩散使得浓度层之间的重叠成为可能，导致形成部分重叠的层状束并且在孔隙尺度上实现均匀化。这些机制控制了体积和表面浓度方差的演变，产生了特征性的代数衰减（$\sigma ^2 \sim t^{-1/2}$），通过数值模拟得到了验证。当Péclet数远小于表征指数混合和剪切引起的变形开始时间尺度的比率时，可以观察到这种动态。在小的、紧凑的固定床反应器中，这种动态是成立的。聚集动态结合了相关的纵向和随机的横向重叠，解释了与纯随机三维混合相比体积方差较慢的衰减。该模型突出了剪切引起的变形在有限时间内的关键作用。尽管在后期指数拉伸在理论上变得相关，但它在这里只起次要作用，因为剪切驱动的机制在对于传输和反应相关的时间尺度上占主导地位。在反应性系统中，流体-催化剂界面的表面浓度均匀化反映了体积混合趋势，表面浓度方差遵循相同的缩放，直到特征反应时间。表面浓度的均匀化遵循Pe缩放$\sigma _s^2 \propto (\gamma ^* Pe)^{-\alpha /2(\alpha +1)}$，这使得剪切引起的混沌重启模型（Aquino等人，参考文献Ben-Noah, Hidalgo和Dentz2023）可以扩展到线性剪切范围之外（$\alpha = 1$）。在混沌轨迹的发散适中的多孔系统中（Lyapunov指数$\lambda \approx 0.1$），Péclet数相当高（$\textit{Pe} = O(10^2)$），并且由于分子扩散（${\textit {Da}}_{\textit {I}}= O(10^1)$），反应相对较慢，强剪切（$\alpha \approx 2$）可以在有限时间内主导混合并控制最终的孔隙尺度表面反应。

致谢
我们感谢欧盟的H2020研究与创新计划，授予编号为951801的资助协议（MOST H2020-EIC-FETPROACT-2019-951801）。计算得到了瑞典国家超级计算学术基础设施（NAISS）提供的资源的支持，部分资金由瑞典研究委员会通过编号为2022-06725的资助协议提供。作者感谢D. Lester就手稿中使用的流体变形计算进行了富有成果的讨论。

利益声明
作者声明没有利益冲突。

附录A
速度梯度张量$\boldsymbol{\epsilon }'(t)$的计算遵循Lester等人（参考文献Lester, Dentz, Le Borgne和de Barros2018）描述的程序。移动的Protean坐标系$\boldsymbol{x}'$与固定的笛卡尔坐标系$\boldsymbol{x}$通过旋转矩阵$\unicode{x1D64C}(t)$相关联，具体关系如下（A1）：
$$
\boldsymbol{x}'=\boldsymbol{x}_0(t)+\unicode{x1D64C}(t)^{{T}} \boldsymbol{x}
$$
其中$\boldsymbol{x}_0(t)$是任意平移向量。变形张量$\unicode{x1D641}(t)$根据速度梯度张量$\boldsymbol{\epsilon }(t) = \{ \boldsymbol{\nabla } \boldsymbol{u}[\boldsymbol{x}(t,\boldsymbol{X})] \}^{{T}}$定义了沿拉格朗日轨迹的变形（A2）：
$$
\frac {\mathrm{d}\unicode{x1D641}}{\mathrm{d} t}= \boldsymbol{\epsilon }(t) \unicode{x1D641}(t)
$$
沿Protean坐标系计算的速度梯度张量为（A3）：
$$
\boldsymbol{\epsilon }'(t)=\unicode{x1D64C}(t)^{{T}}\boldsymbol{\epsilon }(t)\unicode{x1D64C}(t) + \unicode{x1D63C}(t)
$$
该张量由两部分组成，第二部分量化了移动参考系的影响（A4）：
$$
\unicode{x1D63C}(t)= \dfrac {\mathrm{d} \unicode{x1D64C}(t)^{{T}} }{\mathrm{d} t} \ \unicode{x1D64C}(t)
$$
通过应用正确的变换（A1）并利用定义（A3）和（A4），可以将速度梯度张量转换为上三角矩阵，并沿Protean坐标系重新定义（A2）（A5）：
$$
\frac {\mathrm{d}\unicode{x1D641}^{{\kern1pt} \prime}}{\mathrm{d} t}= \boldsymbol{\epsilon }'(t) \unicode{x1D641}^{{\kern1pt} \prime}(t)
$$
通过对$\unicode{x1D641}^{{\kern1pt} \prime}(t)$进行连续积分，可以很容易地求解它，从而明确地确定导致变形的剪切和拉伸分量。Protean坐标系的方向通过旋转矩阵$\unicode{x1D64C}(t)=\unicode{x1D64C}_1(t)\unicode{x1D64C}_2(t)$获得。给定单位旋转轴$\boldsymbol{q}(t)$和角度$\theta (t)$（A6）：
$$
\boldsymbol{q}(t) = \dfrac {1}{\sqrt {v_2^2+v_3^2}} (0, u_3, -u_2)^{{T}}
$$
（A7）
$$
\cos \theta (t) = \dfrac {u_1}{v}
$$
其中$\boldsymbol{v}(t,\boldsymbol{X})= (u_1,u_2,u_3)^{{T}}$是沿计算出的拉格朗日轨迹的瞬时速度，其大小为$v$。第一个旋转矩阵定义为（A8）：
$$
\unicode{x1D64C}_1(t) = \cos \theta (t) \unicode{x1D644} + \sin \theta (t) (\boldsymbol{q}(t)_\times )^{{T}} + (1-\cos \theta (t)) \boldsymbol{q}(t) \otimes \boldsymbol{q}(t)
$$
在上述（A8）中，$\unicode{x1D644}$是单位矩阵，符号$\times$表示叉积，$\otimes$表示张量积。速度梯度张量的第一次旋转产生旋转后的张量$\boldsymbol{\epsilon }^{(1)}(t)$（A9）：
$$
\boldsymbol{\epsilon }^{(1)}(t) =\unicode{x1D64C}_1(t)^{{T}} \boldsymbol{\epsilon }(t) \unicode{x1D64C}_1(t)
$$
第二次旋转通过$\unicode{x1D64C}_2(t)$实现，角度为$\psi (t)$（A10）：
$$
\mathsf{Q}_2(t) = \left ( \begin{array}{c@{\quad}c@{\quad}c} 0 & 0 & 0 \\[3pt] 0 & \cos \psi (t)& -\sin \psi (t) \\[3pt] 0 & \sin \psi (t)& \cos \psi (t) \end{array} \right )
$$
最终的张量由（A3）给出，可以重写为（A11）：
$$
\boldsymbol{\epsilon }'(t) =\unicode{x1D64C}_2(t)^{{T}} \boldsymbol{\epsilon }^{(1)}(t) \unicode{x1D64C}_2(t)
$$
通过选择满足以下常微分方程的$\psi (t)$，可以强制$\boldsymbol{\epsilon }'(t)$满足上三角矩阵的条件（A12）：
$$
\dfrac {\mathrm{d}\psi (t)}{\mathrm{d} t} = a(t) \cos ^2\psi (t) + b(t) \cos ^2\psi (t) +c(t) \cos \psi (t) \sin \psi (t)
$$
其中系数$a(t)$、$b(t)$和$c(t)$的定义如下（A13）：
$$
a(t) = \epsilon ^{(1)}_{32} + \dfrac {u_2}{v+u_1}\epsilon ^{(1)}_{31} - \dfrac {u_3}{v+u_1}\epsilon ^{(1)}_{21}
$$
（A14）
$$
b(t) = \epsilon ^{(1)}_{23} + \dfrac {u_2}{v+u_1}\epsilon ^{(1)}_{31} - \dfrac {u_3}{v+u_1}\epsilon ^{(1)}_{21}
$$
（A15）
$$
c(t) = \epsilon ^{(1)}_{22} + \epsilon ^{(1)}_{33}
$$
（A16）
其中$\epsilon ^{(1)}_{\textit{ij}}$是张量$\boldsymbol{\epsilon }^{(1)}$的$i,j$分量。Lester等人（参考文献Lester, Dentz, Le Borgne和de Barros2018）指出，微分方程（A12）虽然允许多个解，但在初始角度条件$\psi (0)$的特定值下存在一个吸引性轨迹。该值可以通过最小化$\mathrm{d} \dot {\psi } / \mathrm{d} \psi$在短时间内估算出来，其中$\dot {\psi }=\mathrm{d}\psi /\mathrm{d} t$。一旦确定了$\psi _0$，方程（A12）的解就提供了沿拉格朗日轨迹的$\psi (t)$、$\unicode{x1D64C}_2(t)$的时间解。根据（A11），可以依次计算变形张量$\unicode{x1D641}^{{\kern1pt} \prime}(t)$（A16）：
$$
F'_{11} = \dfrac {v(t)}{v(0)}
$$
（A17）
$$
F'_{22} = \exp \bigg [ \int _0^t \epsilon '_{22}(t)\ \mathrm{d} t \bigg ]
$$
（A18）
$$
F'_{33} = \exp \bigg [ \int _0^t \epsilon '_{33}(t)\ \mathrm{d} t \bigg ]
$$
（A19）
$$
F'_{12} = v(t) \int _0^t \dfrac {\epsilon '_{12}(t) F'_{22}(t)}{v(t)}\ \mathrm{d} t
$$
（A20）
$$
F'_{23} = F'_{22}(t) \int _0^t \dfrac {\epsilon '_{23}(t) F'_{33}(t)}{F'_{22}(t)}\ \mathrm{d} t
$$
（A21）
$$
F'_{13} = v(t) \int _0^t \dfrac {\epsilon '_{12}(t) F'_{23}(t)+\epsilon '_{13}(t) F'_{33}(t)}{v(t)}\ \mathrm{d} t
$$
图14：上部面板显示了速度梯度张量对角分量$\epsilon '_{\textit{ii}}$的概率分布函数$P$，对于特征去相关时间$t'_{\textit {num}}=10^3$（a）和$t'_{\textit {num}}=10^4$（b）的VH流动情况。平均值$\langle \epsilon '_{22}\rangle$分别对应Lyapunov指数$\lambda =0.48$（a）和$\lambda =0.57$（b）。下部面板显示了速度自相关函数$\langle v'(t)v'(0)\rangle /\langle v'^2(0)\rangle$（其中$v'(t) = v(t)-\langle v(t) \rangle$表示速度幅值波动），用于识别去相关特征时间$t'_{\textit {num}}=10^3$（a）和$t'_{\textit {num}}=10^4$（b）。红色虚线表示指数衰减$\langle v'(t)v'(0)\rangle /\langle v'^2(0)\rangle = \exp (-t_{\textit {num}}/t'_{\textit {num}})$。右下角面板展示了$R/d=4.7$和$R/d=3.1$的-reactor模拟的自相关函数。右面板（c）显示了在VH流动上计算出的轨迹。我们在Lester等人（参考文献Lester, Dentz, Le Borgne和de Barros2018）描述的变量螺旋度（VH）流动中测试了上述过程。结果如图14所示，右面板（c）展示了一个通过拉格朗日跟踪算法计算的典型轨迹。我们特别关注轨迹求解器的时间分辨率对新结果精度的影响。一旦选择了时间分辨率，VH流动支持的轨迹会在某个特征时间$t'_{\textit {num}}$上发生去相关，见图14的下部面板。我们发现，对于时间分辨率$t'_{\textit {num}} = 10^3$和$t'_{\textit {num}} = 10^4$，对应的平均对角分量分别为$\langle \epsilon '_{\textit{ii}} \rangle = (0, 0.48,-0.48)$和$\langle \epsilon '_{\textit{ii}} \rangle = (0, 0.57,-0.57)$。后一个值与Lyapunov指数的理论最大值$\ln (2)$非常接近，表明$t'_{\textit {num}}= 10^4$的时间分辨率是良好的。相同的分辨率也用于图14下部面板（b）中显示的$R/d=4.7$和$R/d=3.1$的 reactors轨迹模拟。