本文档是荷兰代尔夫特理工大学（Delft University of Technology）的研究人员Cees Haringa、Henk J. Noorman和Robert F. Mudde于2017年在《Chemical Engineering Science》期刊上发表的一篇原创性研究论文。该研究属于化学工程与生物过程工程交叉领域，具体聚焦于计算流体动力学（CFD）与反应动力学（RD）耦合模拟在（生物）化学反应器设计与优化中的应用。

一、 研究团队与发表信息

本项研究的主要作者为Cees Haringa（第一作者及通讯作者）、Henk J. Noorman和Robert F. Mudde。他们分别来自代尔夫特理工大学的化学工程系传递现象课题组（Transport Phenomena）、DSM生物技术中心（DSM Biotechnology Center）以及生物技术系的生物过程工程组（Bioprocess Engineering）。这项研究成果于2016年7月被接收，并于2017年正式发表在化工领域权威期刊《Chemical Engineering Science》上，卷号为157，页码为159-168。

二、 学术背景与研究目标

在大型工业规模的（生物）化学反应器中，由于混合时间（bulk circulation time）可能与甚至长于底物代谢周转时间（turnover time），反应器内会形成显著的空间底物浓度梯度。对于悬浮的微生物或催化剂颗粒而言，这种空间梯度会转化为其经历的时间上的底物浓度波动。这种动态环境会影响微生物的代谢速率（如比生长率μ和比生产速率qp），导致群体异质性，并使得在理想混合的小试规模反应器中获得的微生物性能数据无法准确预测其在非理想混合的工业规模反应器中的表现。为了在工艺和菌种开发阶段就考虑这种梯度效应，通常使用缩小规模模拟器（scale-down simulators）或计算模拟方法进行研究。

其中，将计算流体动力学（CFD）与反应动力学（RD）耦合是一种强大的模拟工具。在耦合方法中，欧拉-拉格朗日（Euler-Lagrange, EL）方法因其独特优势而受到关注。该方法将生物相离散为大量携带内部状态参数（如胞内代谢物浓度）的虚拟粒子（文中称为parcels，包裹），从微生物的视角追踪其在反应器内的运动及所经历的环境变化，从而能够方便地研究复杂的、非平衡的代谢动力学响应。然而，EL方法面临巨大的计算负担挑战。计算成本主要受两个参数影响：1）模拟所需的虚拟粒子数量（Np）；2）模拟的时间步长（Δt）。粒子数量不足会导致生物质浓度在空间分布上出现人为的虚假波动，进而引发虚假的底物浓度梯度，影响模拟精度。而时间步长过大则会影响轨迹计算的准确性并导致质量不平衡。因此，如何以最小的计算成本获得准确可靠的模拟结果，是该方法实际应用的关键瓶颈。

本研究的目标正是为了解决这一问题。作者旨在系统地研究Np和Δt这两个关键参数对EL模拟结果的准确性和计算时间的影响，并基于此建立一套实用的、可预测的指导原则，帮助研究者在进行工业规模生物反应器的EL-CFD模拟时，能够预先确定所需的最小粒子数和最大允许时间步长，从而在保证结果准确性的前提下，最大限度地提高计算效率。

三、 详细研究流程与方法

本研究的工作流程主要包括理论分析、数值模拟实施、参数化研究以及结果验证与讨论几个核心环节。

第一环节：理论分析与预测模型建立。 本研究首先从理论层面分析了由离散粒子分布引起的“人工浓度梯度”问题。研究者将一个具有代表性的虚拟粒子及其周围的影响体积（Vp）抽象为一个球形区域。由于粒子是点状的，底物消耗仅发生在包含该粒子的网格单元内，而周围流体的底物浓度较高，从而在Vp内形成了一个稳态的浓度梯度。他们定义了一个无量纲参数β来量化这个人工梯度的幅度（β = (Cs,b - Cs,p)/Cs,p，其中Cs,b和Cs,p分别为Vp边界和粒子处的底物浓度）。

通过分析混合时间尺度（τm,p，基于湍流扩散理论）与反应时间尺度（τr,c）之间的平衡，并假设反应遵循Monod动力学且在底物浓度远小于半饱和常数Ks时简化为一级反应（此时β最大，为最坏情况），研究者建立了β与过程参数之间的关系。他们通过求解一维扩散方程（使用MATLAB），并结合大量数值实验，最终推导出最大人工梯度βm与（Vc/Vp）^(²⁄₃) * (τm,p / τr,c)成正比的经验关系式（公式3）。

考虑到实际模拟中粒子是随机运动的，每个网格单元内的粒子数np,c服从二项分布。因此，整个计算域内由人工梯度引起的浓度波动的总体效应，可以用一个体积平均的浓度变异系数χ（定义为每个网格内浓度时间序列的标准差与均值的比值，再进行体积平均）来量化。研究者假设χ与每个网格的βm,c和粒子数标准差σp,c的乘积的体积平均值成正比，即 χ ∝ < βm,c σp,c >，比例常数为α。通过联立上述关系式，最终推导出预测所需最小粒子数Np的理论公式（公式6）。该公式将Np与反应器体积（Vt）、网格单元体积（Vc）、湍流扩散系数（Dt,c）、动力学参数（最大比反应速率ks,max、半饱和常数Ks、生物质浓度Cx）以及允许的波动阈值χ联系起来。其中，几何相关的比例常数α需要通过基准模拟来确定。

第二环节：CFD模拟实施与参数设置。 所有CFD模拟均在商业软件ANSYS Fluent 15.07中完成。研究选择了一个假设的工业规模搅拌釜反应器作为几何模型：平底罐，高与直径均为5米，采用一个六直叶涡轮（Rushton Turbine, RT）搅拌，并设置了挡板。为了节省计算资源，模拟了1/6的对称区域。使用了多种网格进行测试，包括51k、166k、341k的六面体网格和一个147k的四面体网格。此外，还使用了一个四斜叶涡轮（Pitch Blade Turbine, PBT）的几何模型进行对比，以检验α常数是否受几何构型影响。

流体流动采用稳态的雷诺平均Navier-Stokes（RANS）方法和标准k-ε湍流模型进行模拟。搅拌桨区域采用多重参考系（MRF）模型处理。底物（葡萄糖）的输运通过求解标量输运方程实现。初始时，仅启用欧拉（体积）反应源项，待浓度场达到稳态后，切换到瞬态求解并引入拉格朗日粒子。

拉格朗日相被设置为质量less粒子（因微生物的斯托克斯数St << 0.01），其运动由当地流体速度加上离散随机游走（Discrete Random Walk, DRW）模型模拟的湍流脉动决定。关键的创新在于反应耦合：研究者通过用户自定义函数（UDFs）实现了粒子耦合的反应模型。具体而言，每个粒子根据其所在网格单元的底物浓度（Cs,c），按照Monod动力学计算其消耗速率（rs,p）。这个消耗速率再作为源项（负值）施加到该网格的欧拉相标量方程中。研究测试了两种耦合方式：方式I（非迭代式）在每个时间步开始时计算源项并保持不变；方式II（迭代式）在每个流场迭代步中更新源项。

物理条件设定为水相，底物为葡萄糖，模拟的微生物是产黄青霉（Penicillium chrysogenum），其动力学参数来自文献。生物质浓度设为10 g/L。大多数模拟在搅拌转速为1.026 rps下进行，也测试了0.75 rps以考察功率输入的影响。时间步长Δt通常设为0.01秒，但也系统性地改变了Δt以研究其影响。

第三环节：参数化研究与数据采集。 研究设计了一系列模拟案例（见表1），系统地改变关键参数：虚拟粒子数量Np（从5k到2000k）、时间步长Δt（从0.001秒到5秒）、网格类型和密度、耦合方法（I或II）、以及搅拌桨类型（RT或PBT）。

在每个模拟中，监测了以下关键数据：整个计算域的平均底物浓度、总底物消耗速率、每个网格单元的瞬时底物浓度（Cs,c）和粒子数（np,c）、以及每个粒子“看到”的底物浓度（Cs,p）。为了消除初始瞬态的影响，在启用拉格朗日反应后，先舍弃相当于15个搅拌桨回转时间的模拟数据，然后持续监测60-300个回转时间的数据用于统计分析。主要分析的指标包括：1）人工梯度波动水平χ；2）拉格朗日模拟与纯欧拉模拟（均匀生物质分布）所得平均浓度的一致性；3）欧拉相与拉格朗日相之间的底物质量平衡闭合误差λ。

第四环节：数据分析与模型验证。 通过后处理（主要在MATLAB中完成），将模拟中观测到的χ值与根据网格和模拟参数计算出的预测值<β*m,c σp,c>（β*m,c是缩放后的βm,c）进行对比，以确定公式6中的比例常数α。同时，比较了不同Np下拉格朗日模拟的平均浓度与对应网格下纯欧拉模拟的参考值，以验证EL方法的准确性。还分析了质量不平衡误差λ以及Δt对χ和λ的影响。

四、 主要研究结果

1. 人工梯度波动（χ）与预测模型的关系： 研究结果证实，模拟中观测到的人工浓度波动χ与理论预测值<β*m,c σp,c>之间存在良好的线性关系（见图4），这与理论假设一致。通过线性拟合，确定了比例常数α的值：对于Rushton涡轮几何构型，α约为0.27；对于PBT几何构型，α约为0.19。这一差异被归因于不同几何构型导致的全局流场模式和粒子运动特性的不同。重要的是，α值被证明与网格类型（六面体/四面体）、网格密度、功率输入（在测试范围内）以及时间步长（在合理范围内）无关，也与两种反应耦合方式无关。在双Rushton涡轮和滑动网格（Sliding Mesh）模拟中也得到了相同的α值（0.27）。这一发现表明，α是一个仅取决于反应器几何构型的常数，这极大地增强了预测公式（公式6）的普适性和实用性。

2. 平均浓度的一致性： 随着粒子数Np增加（即<β*m,c σp,c>减小），拉格朗日模拟得到的体积平均底物浓度逐渐收敛，并在<β*m,c σp,c>低于约0.5时趋于稳定（见图5）。这个稳定值与相同网格下纯欧拉模拟（均匀生物质分布）得到的平均浓度基本一致（见表2），验证了EL方法在足够多粒子数下能够准确再现反应器的宏观行为。然而，研究也发现了一个与DRW湍流模型相关的网格依赖性现象：在较粗的网格中，粒子容易在近壁区域（特别是靠近底物加料的顶部区域）异常滞留，导致局部生物质浓度虚高，进而使整体平均浓度略低于欧拉模拟值。网格加密可以缓解此问题。这提示使用者，在检查模拟准确性时，除了关注Np，还需验证粒子是否在特定区域异常聚集。

3. 相间质量平衡： 质量不平衡误差λ被定义为拉格朗日相消耗的底物与从欧拉相移出的底物之间的相对差异。研究发现，当χ < 0.05且Δt满足特定条件时，λ的绝对值始终低于1%（见图6），表明质量平衡闭合良好。两种反应耦合方式（I和II）在误差方面没有表现出显著差异，但方式I（非迭代）的计算成本更低。

4. 时间步长Δt的影响： 研究引入了无量纲参数θ = (Δt * ks,max * Cx) / Ks，用于将Δt与反应时间尺度相关联。结果显示（见图7），只要θ < 1，时间步长Δt对人工梯度波动水平χ几乎没有影响。然而，Δt对质量不平衡误差λ有显著影响：当θ > 0.1时，λ开始显著增大，这主要是由于在一个时间步内，网格单元中的底物可能被过度消耗（“裁剪”效应）。此外，为确保粒子轨迹计算的稳定性，Δt不应超过搅拌桨回转周期的十分之一（即Δt < 1/(10*ns)）。因此，最大允许Δt应由这两个条件共同决定。

五、 研究结论与价值

本研究成功地为高效、准确地进行欧拉-拉格朗日生物反应器模拟提供了清晰、可操作的设置指南。主要结论和指导原则如下：

1. 确定最小粒子数（Np）：可以使用推导出的公式（6）预先计算所需的最小粒子数，以确保人工浓度波动χ低于可接受的阈值（例如χ = 0.05）。对于典型情况，Np的数量级通常在10^5左右。公式中的关键参数α取决于反应器几何形状，对于Rushton涡轮可取0.27，对于PBT可取0.19。
2. 确定最大时间步长（Δt）：为确保计算稳定性和质量平衡，建议时间步长满足两个条件：Δt < min( 0.1 * Ks / (ks,max * Cx)， 1 / (10 * ns) )。第一个条件基于反应动力学，第二个条件基于流体运动。
3. 优化计算设置：为了减少计算时间，建议：（a）采用非迭代的反应耦合方式（即方式I）；（b）使用多重参考系（MRF）方法模拟搅拌桨；（c）将搅拌桨和挡板的边界条件设置为“界面”（interface）而非“壁面”（wall），以防止粒子滞留；（d）谨慎使用离散随机游走（DRW）湍流模型，因其可能导致粒子在低湍流区异常聚集，需通过网格加密或其他更先进的模型（如连续随机游走CRW）来缓解。

本研究的科学价值在于首次从理论和模拟上系统地量化了EL-CFD模拟中两个关键计算参数（Np和Δt）对结果精度的影响，并建立了可预测的定量关系。这为复杂生物过程的多尺度模拟（耦合详细代谢模型与反应器流体力学）提供了重要的方法论基础，使得原本计算成本极高的工业规模模拟变得更为可行。

其应用价值非常直接：为从事生物反应器模拟和放大研究的工程师与科学家提供了一套实用的“操作手册”。遵循这些指南，研究者可以在项目开始前合理配置计算资源，避免因参数设置不当而导致模拟结果不可靠或计算时间过长，从而更有效地利用EL-CFD工具来研究梯度效应、优化反应器设计、指导缩小规模实验的构建，并最终助力于更稳健的工业生物过程开发。

六、 研究亮点

1. 问题导向的实用性研究：本研究并非开发新算法，而是针对EL-CFD这一已有方法在实际应用中的核心痛点——计算成本与精度权衡——进行了深入剖析，成果具有很高的工程实用价值。
2. 理论推导与数值验证相结合：研究从基本的物理和数学原理出发，推导出预测人工梯度的理论框架，并通过大量系统的CFD模拟数据对理论预测进行了验证和校准，使得最终给出的指南有坚实的理论基础和实证支持。
3. 参数影响的系统量化：首次明确给出了Np和Δt的定量选择依据。特别是推导出的Np预测公式，将粒子数与网格属性、流体力学参数和反应动力学参数直接联系起来，是一个重要的创新点。
4. 揭示并警示了模型局限性：研究不仅给出了优化指南，还明确指出了当前常用模型（如DRW湍流模型）可能引入的误差（粒子在壁面聚集），并提出了解决方案，这对模拟实践者具有重要的警示和指导意义。
5. 跨几何构型的参数标定：研究发现并确认了比例常数α的几何依赖性，虽然增加了应用时的一次性标定工作，但使指南更具一般性，不局限于特定搅拌型式。

七、 其他有价值内容

论文附录提供了宝贵的计算时间基准测试数据（附录A）和流场验证结果（附录B）。附录A详细分析了不同模拟设置（仅跟踪、跟踪+反应、跟踪+反应+数据输出）对单粒子每时间步计算时间（tcomp）的影响，指出文件输出是主要耗时环节，并给出了并行计算效率的参考。附录B验证了所使用的RANS k-ε模型和MRF方法在预测所研究搅拌槽的湍流动能、耗散率和功率准数方面与文献实验数据基本吻合，确保了流体力学背景的可靠性，尽管也指出了四面体网格在预测湍流参数方面的不足。这些信息对于实际应用中的计算资源规划和模型可信度评估非常有帮助。