CFD 应用于各个领域，包括生物技术，用于研究、开发和教育。最广泛使用的方法是离散纳维-斯托克斯（Navier-Stokes）方程的有限体积法。此外，这种方法用于大多数商业和开源 CFD 代码（例如 Ansys Fluent 和 CFX、Siemens Simcenter Star-CCM+、OpenFOAM），有限元方法 (Autodesk CFD) 或 Lattice-Boltzmann 方法 (M -Star CFD、OpenLB）也有被使用。在后者中，玻尔兹曼（Boltzmann）方程被离散化，而不是纳维-斯托克斯方程。本章仅讨论有限体积法。为了描述生物反应器中的单相层流，质量和动量守恒方程被离散化并迭代求解。单相考虑仅适用于可以假设恒定流体表面（无涡流形成）的搅拌生物反应器。

对于大多数生物技术应用，需要产生湍流以确保细胞悬浮、氧气分散和营养物质均匀分布的充分混合。然而，对于用于培养可能对流体动力应力敏感的细胞的生物反应器来说，情况并非如此。有多种方法可以解释湍流，其中雷诺（Reynolds）平均纳维-斯托克斯 (RANS) 方法由于其经济优势而被最广泛使用。该方法使用附加模型完整地描述了湍流，并且仅计算平均通量。有多种不同的 RANS 湍流模型。然而，基于涡流-粘性的模型，例如κ-ε-模型 （方程（4）和（5））和κ-ω模型（方程可以在 Wilcox 的文章中找到），是最广泛使用的。如果需要更详细地考虑湍流现象，可以使用非稳态RANS、混合RANS大涡模拟（LES）、LES或直接数值模拟（DNS），但计算复杂度会增加。事实上，DNS 的计算量非常大，以至于它仅用于研究中的简单几何结构。Wilcox 和 Rodriguez 描述了湍流模型的详细概述。

机械驱动生物反应器根据将动力引入系统的方法进行表征，如图 1 所示。对于搅拌系统，主要使用两种不同的方法来模拟搅拌桨的运动。使用冻结转子方法（也称为多参考系方法），搅拌桨不旋转，但在搅拌桨周围的区域中，守恒方程在旋转参考系中求解。因此，不存在计算网格的位移或重新计算，这使得该方法资源高效。然而，该方法仅适用于搅拌生物反应器的稳态考虑。另一种方法是滑动网格方法，其中搅拌桨和搅拌区的网格相对于生物反应器的其余部分移动。该方法的计算强度明显高于冻结转子方法，但允许瞬态方法。轨道振摇和波浪混合生物反应器也有两种不同的方法，它们都使用瞬态方法。计算网格可以通过旋转和/或平移来模拟运动。一种不太常见的方法是仅操纵生物反应器内作用力的方向。例如，Zhu等人利用离心力的操纵成功地模拟了轨道振摇系统。

就工艺工程表征而言，单相模拟可用于确定流场以及非充气单位功率输入，功率输入的确定类似于等式：(1)。可以使用方程式（6）来确定扭矩，而不是测量扭矩。其中Fp和Fτ对应于压力和粘性力，并且τ 为从网格单元中心到旋转轴的距离（方程（7）和（8））。

法向梯度和剪切梯度可以根据速度场或其梯度来确定。Wollny 中描述了详细的推导过程。通过使用RANS湍流模型，湍流能量耗散率ε可以确定。这也可以直接计算，例如，当使用κ-ε--模型时，或者可以从κ和ω确定，即当使用κ-ω-模型时（方程（9））。能量耗散率可用于预测柯尔莫哥洛夫长度λκ（方程（10））。如果对每个网格单元都确定了这一点，则可以确定生物反应器中的柯尔莫哥洛夫长度分布。

搅拌生物反应器中的混合时间也可以根据静态单相模拟来确定。为此，速度场的收敛稳态解被用作基础。随后，虚拟跟踪器Tm 可以通过标量传输方程（方程（11））输入系统。

随后的瞬态观察显示示踪剂如何通过对流和扩散（包括湍流扩散Dt）随着时间的推移在系统中分布。与脱色方法类似，可以用于确定混合时间。

如果正在研究轨道振摇、波浪混合或涡流形成的搅拌生物反应器系统中的流体流动，则除了液相之外还必须对气相进行建模。研究多相系统存在不同的方法。然而，流体体积（VOF）方法适用于连续两相系统，其中可以假设连续性。连续性方程对应于方程 (12)， 而动量方程对应方程（13）。VOF 方法使用混合流体，其中物理特性χ，比如密度ρ和粘度μ按相分数加权（方程（14）和（15））。表面张力σ在这里可以通过使用连续表面力（CSF）模型来考虑。

与单相模拟类似，功率输入也可以使用扭矩值来确定。这里，必须考虑相分数以及这是瞬态确定的事实，可能必须随时间对平均单位功率输入进行平均。剪切梯度和法向梯度以及柯尔莫哥洛夫长度也可以通过与单相模拟类似的方式确定。

这些采用单纯表面通气的系统，KLa也可以确定。与实验研究相比，KL和a在 CFD 模拟中单独确定。对于表面通气系统，可以根据液体表面积与工作体积的商来确定特定界面，其中液体表面积使用等值面确定αwater=0.5。KL值可以通过多种模型来近似。Higbie 的模型被广泛使用（方程（16）），并且需要确定能量耗散率。混合时间可以由虚拟示踪剂确定。然而，在 VOF 模型中，确保虚拟示踪剂不会迁移到气相中非常重要。商业和开源 CFD 代码中使用了不同的模型方法（OpenFOAM 中的相标量传输、Ansys Fluent 中的用户定义标量）。

然而，如果考虑分散系统，VOF 模型就无法表示这种复杂的现象。分散系统包括主动通气系统和带有颗粒（例如微载体）的系统。与 VOF 模型相反，Euler-Euler 模型不使用混合流体，而是单独求解每个相的守恒方程。然而，对于所有相仅计算一个压力场。使用Euler-Euler模型，还可以考虑界面力，例如阻力、升力和虚拟质量力，这使得可以对定义尺寸的颗粒进行建模（不同的颗粒尺寸可以通过附加相进行建模， 但这是一种计算密集型方法）。Seidel 、Pourtousi和Chuang及Hibiki等人详细描述了这些力及其计算。如果已知气泡直径d可以假设恒定，单位界面a也可以用 Euler-Euler 模型确定（方程（17），对于单分散系统：d32=d）。

Seidel 和 Eibl 已经能够证明，在 3 L 搅拌生物反应器中，在相同的模拟条件下，Kla值可以描述为气泡直径的负指数幂函数。然而，将 CFD 与群体平衡建模 (PBM) 相结合可以最大限度地减少模拟Kla值对所选气泡直径的严重依赖性。这是通过对具有气泡尺寸分布的气泡群进行建模来实现的。此外，PBM还可以通过描述气泡直径的数密度函数的变化来模拟气泡破裂和聚结等复杂现象，源项由四个部分组成：聚结消失、聚结产生、破裂消失，破裂产生。种群平衡方程闭合模型的详细描述和介绍可以在 Liao 和 Lucas的文章中找到，Seidel 和 Eibl的文章还描述了聚结和破裂模型的选择如何影响KLa值。最常用的方法包括 OpenFOAM 中使用的类方法和 Ansys Fluent 中使用的矩量方法。潜在方法的详细概述可以在 Nguyen 和Li等人的文章中找到，Seidel 等人则讨论了它们的优点和缺点。

Euler-Euler 模型还可用于执行悬浮研究并确定NS1标准。然而，关于该标准的 CFD 研究还不是很广泛。为了将微载体建模为固相，需要颗粒流模型（Granular Flow Model，GFM），颗粒流动力学理论（KTGF）方法被广泛用于此目的。出于建模目的，引入了颗粒压力，这需要颗粒相动量方程的梯度项。为了计算颗粒压力，还需要计算颗粒温度，它影响颗粒的碰撞行为。Syamlal 等人的文章中可以找到如何对颗粒压力进行建模的示例。使用 Euler-Euler 模型时，不是对单个颗粒（气泡或微载体）进行建模，而是仅对它们的相分数进行建模。然而，如果需要对单个颗粒进行建模，则应使用Euler-Lagrange模型，该模型使用常微分方程描述颗粒的路径。然而，Euler-Lagrange模型没有考虑湍流对颗粒路径的影响，这可以使用离散随机游走模型来解释。事实上，与 Euler-Euler 模型相比，Euler-Lagrange 模型在生物技术应用中的使用频率较低。其主要原因是计算时间随着颗粒数量的增加而线性增加。Zieringer 和 Takors 建议使用Euler-Lagrange 模型最大颗粒馏分为10%。因此，与Jossen的研究一样，使用10g/L (1% 颗粒馏分）微载体，大约需要对40x10^6 载体/L单独建模。

根据所研究的生物反应器系统和工艺参数，应考虑不同的模型，以避免不经济的 CFD 模拟。应该指出的是，无论 CFD 模拟使用什么模型，它都永远无法真实反映现实。除了建模错误之外，离散化错误、迭代收敛错误、舍入错误、编程错误和用户错误都可能导致不准确。舍入误差通常可以忽略不计（双精度浮点数），而用户错误理想情况下应该完全避免。离散化误差源自网格单元的空间离散化和时间步长的时间离散化。用于估计因网格选择而产生的离散化误差的一种广泛使用的方法是Richardson外推法或相关的网格收敛指数，它使用不同分辨率的网格，考虑了目标参数的变化（例如单位功率输入或混合时间）。Ramírez和 Seidel 等人描述了此方法的详细应用。或者，也可以使用网格系统细化或曲线拟合方法等策略。

如果估计的离散化误差已知，则可以选择经济上可接受的网格，因为随着网格单元数量的增加，计算时间呈指数增加，而离散化误差渐近地接近零。对于瞬态模拟，空间和时间离散都发挥着作用。如果选择的时间步长太长，可能会导致模拟不稳定和不准确。与更精细的计算网格一样，随着时间步长变短，计算时间也会线性增加。Courant-Friedrichs-Lewy (CFL) 数有助于估计时间步长。如果使用显式时间积分方案，高于 1 的 CFL 值可能会导致非物理值和发散。即使在隐式方案中，CFL 数也最好低于 1，以避免出现大错误。CFL 数是针对每个网格单元单独计算的，在单个生物反应器模型中可能会有很大差异。因此，任何单个网格单元的最大 CFL 数量不应超过 1。为了确定CFD结果的总误差，需要通过物理实验进行验证。

验证

CFD 模拟的验证非常重要，美国航空航天学会将其定义为“从模型的预期用途的角度确定模型准确表示现实世界程度的过程”。为了进行验证，可以进行不同的实验研究。应该指出的是，实验也并非没有测量误差，统计报表中也必须考虑到这一点。根据目标过程变量，不同的方法适合验证不同的生物反应器模拟。最基本的因素之一是流动模式。颗粒图像测速 (PIV)、激光多普勒流速测定和激光诱导荧光适用于测量生物反应器中的速度，特别适用于搅拌生物反应器。对于轨道振摇和波浪式混合系统，通常仅测量液体高度或分布。

如果功率输入相关，则将其用作验证标准是合理的。因此，围绕旋转轴的扭矩测量适用于搅拌、轨道振摇和波浪式混合生物反应器。或者，可以通过测量提供给电机的电流或通过量热法来确定单位功率输入。Villiger等人描述了一种方法，该方法允许使用聚（甲基丙烯酸甲酯）聚体来确定最大流体动力应力，该方法也适用于验证目的。此外，混合时间也可以作为验证参数。基于 CFD 模拟中是全局还是局部评估混合时间，根据 DECHEMA e.V. 一次性技术工作组，脱色方法或电导率测量都适合验证。然而，这些并不是唯一有效的方法，例如，Vivek 等人使用拉曼光谱探头通过实验确定混合时间，以验证 CFD 确定的混合时间。

对于通气系统，要么Kla值可以直接用作验证参数，也可以使用气泡尺寸分布。Kla可以使用 Seidel 和 Eibl 或 Garcia-Ochoa 和 Gomez 文章中描述的不同方法来测量，气泡尺寸分布可以通过阴影摄影术、光学多模在线探针、毛细管抽吸探针技术、SOPAT 内窥镜检查或聚焦光束反射率检测来测量。气含率通常由压差法或雷达液位计光学测定，也适用于通气系统。

对于微载体或其它颗粒的悬浮研究，建议使用悬浮行为或NS1验证标准。为此，将要悬浮的颗粒添加到未搅拌的生物反应器中。然后速度逐步增加，直到NS1达到标准。评估要么通过眼睛执行，要么用相机记录。Delafosse和 Loubière 等人通过使用安装在生物反应器一侧的光源和另一侧的相机测量光衰减来改进了这种评估方法。大多数作者都记录了侧面图像以及从下观察的图像，后者使用镜子。为了验证他们的 CFD 模拟，Zhang 等人用高速摄像机记录的图像。然而，他们没有检查在微载体上生长的细胞，而是检查固定化的乳杆菌细胞。

硬件和软件

为了获得最终结果，需要几个步骤，可以分为预处理、处理和后处理。图 2 更详细地解释了各个步骤。通常，当今的 CFD 软件解决方案中包含其中的几个步骤。使用 Ansys 软件解决方案，可以创建几何体 (SpaceClaim)、生成网格 (Meshing)、执行计算 (Fluent)，并使用 CFD-Post 完成后处理。然而，通常使用单独的软件来生成几何图形，例如 Autodesk Inventor、SolidWorks、Solid Edge 或 Salome，以及用于后处理步骤，例如 Tecplot 和 Paraview。最广泛使用的生物反应器建模 CFD 软件解决方案是 Ansys Fluent 和 CFX（商业）以及 OpenFOAM（开源）。然而，西门子的 Simcenter Star-CCM+、Autodesk CFD、COMSOL Multiphysics 和 M-Star CFD 也被使用，后者被宣传为专用于生物反应器应用（所有商业软件）。

图 2. 从问题定义到验证的流程步骤的可视化。以 Infors AG 的 Minifors 2 搅拌生物反应器为例。

除了软件之外，硬件对于模拟的经济性也至关重要，硬件性能、购买价格和功耗都发挥着作用。当前所有 CFD 软件解决方案都允许并行计算。为此，计算网格被划分为不同的域。然后，各个处理器在各个域中执行计算。像并行计算中通常那样，通过消息传递接口（MPI）或其它接口来调节域之间的通信。创建分区存在不同的算法，其自动化程度和分区时间不同。例如，Pellegrini 和 Roman的 Scotch 算法基于对偶递归双分区并在 OpenFOAM 中实现，只需要域的数量来执行该过程。然而，该算法比自动化程度较低的算法需要更长的分区时间。

虽然并行化减少了所需的计算时间，但更高程度的并行化也需要在处理器之间交换更多数据，这对计算时间产生负面影响。图 3 显示了 3 L 搅拌生物反应器建模的相对模拟时间如何随着并行度的增加而变化（interFOAM、OpenFOAM v9）。利用 Seidel 和 Eibl 文章中描述的硬件设置，相对模拟时间可以描述为幂函数，其指数为-0.817。如果没有因通信造成的损失，则预期指数为-1.。Harasek等人还进行了多达 1024 个核的并行化研究，其中可以确定的指数为-0.93。根据Haddadi等人的说法，并行化应该以每个域有 50,000 到 100,000 个单元的方式进行，随着模型复杂性的增加，单元的数量应该减少。

图 3. 相对模拟时间取决于所使用的处理器。使用 OpenFOAM v9 检查的搅拌生物反应器作为测试用例（VOF 模型，瞬态观察 10 s）。所有模拟均进行三次以捕获时间方差。在每种情况下，均使用 Scotch 算法进行分解。Seidel 和 Eibl 更详细地描述了硬件设置。

在表 1 列出的文章中，只有三分之一的作者对所使用的硬件做出了声明。计算中使用了 4到 504 个核。只有五位作者表示他们使用过 HPC 系统，其余作者使用的是台式机。Ansys Fluent和Siemens Simcenter Star-CCM+的当前版本强调使用图形处理单元（GPU）而不是中央处理单元进行计算。这里也可以使用多个 GPU。Ansys 和Siemens的基准测试表明，使用 GPU 可以缩短模拟时间，同时降低采购成本和功耗。在表 1 列出的作者中，只有使用 M-Star CFD 的作者表示使用了 GPU（M-Star CFD 完全依赖于 GPU）。默认情况下，COMSOL、Autodesk CFD、Ansys CFX 和 OpenFOAM 无法使用 GPU 进行计算。然而，有一些修改版本（例如 MixIT）允许 OpenFOAM 使用 GPU。目前的发展表明，经典 CFD 模拟可以被其它技术补充或取代。机器学习技术可用于加速模拟或改进湍流建模，并且物理信息神经网络 (PINN) 的使用越来越多，因为它们能够以 200 倍的速度以相同的精度执行计算。另一种技术是量子 CFD (QCFD)，目前尚未真正应用于生物反应器建模。在量子计算中，系统可以同时处于多种状态的叠加。通过将传统算法应用于量子计算，可以实现速度的巨大提高，并且湍流等模型的使用将变得过时。

表1. 使用CFD的生物反应器工艺表征文献总结。

按功率输入将生物反应器分为搅拌式、波浪混合式和轨道振摇式。3BEE：3叶象耳搅拌桨；3BSS：3叶片分段搅拌桨；A310：水翼A310；MI：船用搅拌桨；P：桨式搅拌桨；PB：斜桨搅拌桨；RT：Rushion涡轮搅拌桨；SD，特殊设计；ST，Smith涡轮搅拌桨；标签后的数字表示有多少搅拌桨叶。

总结和展望

CFD是一种广泛使用的工具，多年来一直在生物技术中发挥作用，并不断得到普及。它可以用来确定功率输入、混合时间和氧传质等经典参数，也可以利用空间和时间分辨率进行进一步的评估。CFD也是确定流体动力应力的一种很好的方法。然而，根据生物反应器系统和参数的不同，应该选择不同的建模方法来确定这些参数。调查显示，在生物反应器建模方面，Ansys公司的Fluent商业软件占主导地位，其次是开源解决方案OpenFOAM。很少有作者在GPU上计算，但是他们仍然使用CPU。即使与最新的软件解决方案结合使用，GPU的成本更低，功耗更低，计算速度更快。尽管目前有相当大的计算能力，但在与工业相关的模拟(如生物反应器)中，湍流和气泡破裂等现象仍然使用半经验模型进行近似模拟。未来，用于QCFD的更强大的GPU甚至量子处理器将提高计算能力，而像PINN这样的方法将加快计算速度，有可能消除对这些模型的需求，并执行更准确、更快的模拟。

原文：S.Seidel, C.Schirmer, R.W.Maschke, et al., Computational Fluid Dynamics for Advanced Characterisation of Bioreactors Used in the Biopharmaceutical Industry: Part I: Literature Review. Computational Fluid Dynamics - Recent Advances, New Perspectives and Applications, 2023.

撰稿人 |  生物工艺与技术

责任编辑 | 胡静

审核人 | 何发