一种基于浸入边界法的搅拌反应釜模拟方法

摘要

本发明提供了一种基于浸入边界法的搅拌反应釜模拟方法，所述方法为：准备阶段，包括计算参数的测定和设置、生成流体网格文件、生成拉格朗日标记点信息和准备计算资源；数值计算，设定计算区域的边界条件并基于流体网格求解速度场与温度场；结合离散单元法可以实现颗粒流体两相流的界面解析模拟，颗粒表面处直接采用浸入边界法施加无滑移边界条件；后处理阶段，模拟完成后输出每个欧拉网格点的信息以及每个颗粒的信息。该方法依据不同边界的特点，在计算中采用混合浸入边界法进行处理并结合大涡模拟与高性能并行计算技术实现搅拌反应釜内的湍流流动大规模数值模拟。该模拟方法可以应用于搅拌反应釜内湍流流动模拟与固体颗粒悬浮的模拟等。

说明书

技术领域

本发明属于搅拌反应釜的数值模拟领域，涉及一种基于浸入边界法的搅拌反应釜模拟方 法，尤其涉及一种基于浸入边界法的实现搅拌反应釜内湍流流动模拟以及搅拌反应釜内固体 悬浮颗粒的模拟的方法。

背景技术

在工业应用中经常需要进行流体的搅拌与混合。搅拌反应釜常应用于矿业、冶金、生物、 石油、食品、制药等化工领域。由于搅拌反应釜中存在旋转的搅拌桨、复杂的搅拌釜壁面， 甚至还有换热蛇管以及其他的附件等。搅拌反应釜模拟主要的困难在于搅拌桨与反应釜壁面 的相对运动以及搅拌桨复杂几何形状的处理。当前商业软件中常采用的方法主要为多重参考 系方法与滑移网格方法，这两种方法均为贴体网格方法。

Luo等人(PredictionofImpellerInducedFlowsinMixingVesselsUsingMultipleFramesof Reference，EighthEuropeanConferenceonMixing，1994，(136)：549-556.)对内外迭代方法 进行了简化，提出了多重参考系模型(MultipleReferenceFrame，MRF)。在该模型中内外部 区域不再重叠，计算的结果在内外部界面处直接进行匹配。该方法中，内外部界面的位置不 是随意的，需要选在一个流动参数随着时间与周向角度变化较小的位置处。相较于内外迭代 法，多重参考系方法的计算量大大降低。正是因为这个原因，多重参考系模型已经嵌入到大 多数的商用CFD软件中。由于该方法是一个稳态近似，当边界上的流动区域几乎均匀混合或 者搅拌桨与挡板的交互作用较弱时可以使用多重参考系模型。但是，如果需要精确的模拟强 烈的转子和定子间相互作用或者是当搅拌桨与其他静止附件有相互嵌入时，则该方法是不适 用的。

在一般情况下，搅拌桨与挡板的相互作用比较强烈，搅拌槽的湍流特征比较明显，如果 采用稳态模型不能够精确捕捉到这些流动特征。为此，Luo等人(FullFlow-FieldComputation ofMixinginBaffledStirredVessels.ChemEngResDes，1993，71(A3)：342-344.)提出了滑 移网格方法(SlidingMesh，SM)，该方法同多重参考系方法类似，需要划分成静止区域以及 运动区域。但不同于内外迭代法与多重参考系模型，在滑移网格方法中，运动区域内的流动 方程仍然是在惯性坐标系下求解，只是该区域的网格随着搅拌桨的转动而转动。这种方法的 优点是能够较为准确的捕捉搅拌槽内的瞬态现象，并且边界条件设置相对容易。但是该方法 计算量巨大，另外为了考虑内部区域网格的旋转而需要显式的在动量方程中添加加速度项， 并且在界面处的插值运算会严重影响计算效率以及计算精度。

在传统的MRF或者SM方法中，常采用非结构化贴体网格，当搅拌釜以及搅拌桨几何结 构较为简单时网格还比较容易生成。但是当搅拌桨结构复杂或者有其他内构件时(如换热蛇 管等)，生成质量较好的贴体网格非常耗时且非常困难。在有些情况下可能要验证不同桨型的 效果，有些情况下则需要重新划分网格。所述MRF或者SM方法计算区域被划分成静止区域 与运动区域，不同区域间需作差值运算，大大降低了计算精度。当搅拌桨与其他静止附件有 相互嵌入时，这些方法不再适用。

浸入边界法(ImmersedBoundaryMethod，IBM)是通过在控制方程中添加源项的方式来 考虑运动边界对流体的影响，并且对固体边界的运动没有任何的限制，不像滑移网格方法那 样需要一个固定区域与转动区域的界面，而且流体的求解依然在固定的正交网格上进行。因 而，近年来IBM被用于搅拌反应釜的模拟。

Eggels(DirectandLarge-EddySimulationofTurbulentFluidFlowUsingthe Lattice-BoltzmannScheme.IntJHeatFluidFl，1996，17(3)：307-323.)首先实现了在固定正 交网格上进行搅拌反应釜的瞬态模拟。在他的工作中，基于格子Boltzmann方法(Lattice BoltzmannMethod，LBM)来研究流体流动，搅拌桨的运动则通过在动量方程中添加体积力 的方式来实现。模拟得到了精确的平均速度以及湍流脉动量，展现出IBM用于搅拌槽模拟的 可行性以及巨大的潜力。

Verzicco等人(FlowinanImpeller-StirredTankUsinganImmersed-BoundaryMethod.Aiche J，2004，50(6)：1109-1118.)采用IBM对无挡板搅拌槽在中等雷诺数Re下进行了直接数 值模拟(DirectNumericalSimulation，DNS)。由于不存在挡板，所以为了提高计算效率，流 体的求解在柱坐标系下进行。作者同时还指出IBM对于化工领域中更为复杂的问题具有更大 的优势，比如在处理挡板与搅拌桨的相对运动时，通常采用的滑移网格方法需要网格的变形 或者移动，计算量较大，而采用IBM时的计算量与搅拌桨静止情形下大致相同。在该方法中， DNS能够分辨湍流各个尺度的特征，但其总的计算量与Re³成正比，对于当前的计算机条件 还很难实现较高Re的湍流直接数值模拟。

Derksen与VandenAkker(LargeEddySimulationsontheFlowDrivenbyaRushtonTurbine. AicheJ，1999，45(2)：209-221.)在格子Boltzmann方法框架下采用IBM对Rushton桨进 行了大涡模拟。搅拌桨以及釜壁均采用基于Goldstein等人(ModelingaNo-SlipFlowBoundary withanExternalForce-Field.JComputPhys，1993，105(2)：354-366.)提出的直接力方法来 施加边界条件。作者同时指出采用IBM处理复杂运动边界具有很大的灵活性，比如当需要验 证一个新型的搅拌桨或者反应器设计时，只需要生成一套几何表面的标记点，而不需要重新 划分计算网格。但是其仅仅采用一种浸入边界法处理所有的复杂边界问题，没有充分发挥出 浸入边界法的潜力。

Derksen(HighlyResolvedSimulationsofSolidsSuspensioninaSmallMixingTank.AicheJ， 2012，58(10)：3266-3278.)研究了搅拌反应釜内的固体颗粒悬浮，所有复杂壁面包括颗粒 的表面均采用IBM施加无滑移边界条件。这应该是首次实现搅拌反应釜内颗粒悬浮的全尺度 模拟，但其所采用的搅拌槽十分简单且尺寸较小，并且缺乏实验数据的验证。

上述方法中不论是IBM，还是基于IBM的大涡模拟，仅仅是采用一种浸入边界法处理所 有的复杂边界问题，而且很难实现较高Re的湍流直接数值模拟，并且没有一个完整的适用于 搅拌反应釜内湍流流动模拟、搅拌反应釜内固体悬浮颗粒模拟以及搅拌反应釜大规模模拟的 方案，限制了其在工业中的应用。

发明内容

针对上述现有技术中存在的不足，本发明提供了一个适用于搅拌反应釜内湍流流动模拟、 搅拌反应釜内固体悬浮颗粒模拟以及搅拌反应釜大规模模拟的搅拌反应釜模拟方法。该方法 依据不同边界的特点，在计算中分别采用不同的浸入边界法(混合浸入边界法)进行处理并 结合大涡模拟与高性能并行计算技术实现搅拌反应釜内的湍流流动大规模数值模拟。

本发明提供了一种基于浸入边界法的搅拌反应釜模拟方法，所述方法包括以下步骤：

(1)准备阶段，包括计算参数的测定和设置、生成流体网格文件、生成拉格朗日标记点 信息和准备计算资源；

(2)数值计算，设定计算区域的边界条件并基于流体网格求解速度场与温度场，其中计 算区域中运动或者静止边界处的流固耦合采用浸入边界法实现；

(3)结合离散单元法实现颗粒流体两相流的界面解析模拟，颗粒表面处直接采用浸入边 界法施加无滑移边界条件；

(4)后处理阶段，模拟完成后输出每个欧拉网格点的信息以及每个颗粒的信息。

其中，步骤(3)是根据实际操作的需要，选择将数值计算以及其与离散单元法的结合应 用于搅拌反应釜内湍流流动模拟或搅拌反应釜内固体悬浮颗粒的模拟中。

步骤(4)的后处理可以采用Tecplot、Paraview以及Origin等工具进行，同时也可以写 一些小的程序用于辅助数据的处理。

优选地，步骤(1)中计算参数包括数值参数和物理参数。

优选地，所述数值参数包括时间步长、网格数和进程数。

其中，时间步长需要满足一定的条件才能保证数值计算的收敛，在计算中时间步长采用 如下公式进行确定：时间步长＝0.1×网格宽度/(搅拌桨半径×搅拌桨旋转角速度)。

网格数需要根据不同的计算条件进行设置，需要考虑的因素有颗粒的大小、湍流的影响 以及计算量的大小等。在计算中，网格的网格宽度通常为颗粒直径的1/16以保证颗粒表面的 分辨率；另外需要根据雷诺数的大小以及计算量的大小确定合适的网格宽度，当网格宽度不 足以直接模拟湍流时，需要采用湍流模拟，如大涡模拟等。

进程数是根据网格数进行划分的，为了保证计算速度，通常在一个进程中的网格数为128 ×128×128，每个方向的进程数大致为N/128，其中N为其中一个方向上的网格数。

优选地，所述物理参数包括搅拌釜结构参数、流体的密度、流体的粘度、颗粒直径和颗 粒密度。

优选地，搅拌釜结构参数包括搅拌釜的高度、搅拌釜的直径、搅拌桨的几何形状和搅拌 桨的尺寸。

优选地，搅拌釜结构参数采用长度测量工具测定；准确的搅拌釜结构参数可以从其设计 图纸中获得，在没有图纸的情况下可以采用长度测量工具测定，如千分尺等进行长度的测量。

优选地，流体的密度采用密度计进行测量。

优选地，流体的粘度采用毛细管粘度计、落球黏度计或同轴圆筒粘度计中任意一种进行 测量，其根据需要选择合适的工具进行测量。

优选地，颗粒直径采用激光粒度仪进行测量。

优选地，颗粒密度采用密度仪进行测量；所述密度仪主要是粉末和微小颗粒领域的密度 仪。

优选地，步骤(1)中流体网格文件包含每个网格的位置点坐标、速度的初始数值、压力 的初始数值和网格固含率数据。

优选地，所述流体网格为欧拉网格，其可以采用很简单的小程序生成。

优选地，步骤(1)中拉格朗日标记点信息包括搅拌桨表面的标记点信息。

优选地，步骤(1)中拉格朗日标记点信息还包括颗粒表面标记点的信息。在进行搅拌反 应釜内固体悬浮颗粒的模拟时，需要给出颗粒表面标记点的信息。需要指出的是，在搅拌反 应釜模拟过程中需要尽量保证拉格朗日标记点之间的间距与欧拉网格的宽度大致相等。

优选地，步骤(1)中生成拉格朗日标记点信息采用权利要求3中所述的生成方法进行生 成。

优选地，步骤(1)中设定计算程序参数包括数值参数和物理参数。

优选地，所述数值参数包括时间步长、模型系数和进程数。

优选地，所述物理参数包括密度、质量和粘度。

由于模拟过程中，计算区域为矩形区域，故在设定程序参数过程中可根据具体情况将矩 形的六个面设置为无滑移边界，或将顶部平面设置为滑移壁面。

优选地，步骤(1)中准备计算资源包括CPU资源和GPU资源。以确保在并行计算中有 足够的CPU以及GPU计算机资源。

优选地，步骤(2)数值计算阶段包括以下步骤：

(a)读取准备阶段生成的流体网格文件和计算参数，并将其由CPU端传递到GPU端；

(b)在不考虑控制方程中体积力源项的情况下求解纳维-斯托克斯方程得到预估速度；

(c)求解压力泊松方程得到压力修正值，进一步更新预估速度；

(d)流固耦合处理，采用直接力方法处理运动边界得到作用力f₁，采用体积分数方法处 理静止边界得到流固界面处的体积力f₂，将f₁与f₂累加得到作用力f并更新流场速度；

(e)通过作用力f和流场速度更新浸入边界的位置与速度；

其中，纳维-斯托克斯方程是描述粘性不可压缩流体动量守恒的运动方程，简称N-S方程。

流固耦合处理阶段，需要考虑搅拌桨、反应釜壁面以及内部颗粒等对流场的影响。对运 动边界(如运动的搅拌桨与颗粒等)采用直接力方法进行处理，对于静止边界(如搅拌釜壁 面)采用体积分数方法进行处理。

优选地，步骤(d)中采用直接力方法处理运动边界包括以下步骤：

首先，通过狄拉克函数把欧拉网格上的速度u转换到拉格朗日标记点处的U；然后计算 拉格朗日标记点处U的作用力；再将拉格朗日标记处的作用力散布到欧拉网格上得到作用力 f₁；

其中，拉格朗日标记点处的U为拉格朗日标记点处U 的作用力是通过预期速度V^d与当前计算速度的差值计算得到，即其 中RHS包含了对流项、压力项以及粘性项；将拉格朗日标记处的作用力散布到欧拉正交网格 上得到作用力f₁，$f_1\left(x\right)=\underset{l=1}{\overset{N_l}{\Sigma }}F\left(X_l\right){\delta }_h(x-X_l){\mathrm{\Delta V}}_l.$

优选地，步骤(d)中采用体积分数方法处理静止边界是通过网格固含率插值得到流固界 面处的体积力f₂；其中，f₂＝α(u_d-u)/Δt。

步骤(e)中更新浸入边界的位置与速度包括以下三种情况中任意一种或至少两种：

情况①：搅拌反应釜壁面始终保持静止，不更新网格固含率，此部分不需要作任何处理；

情况②：更新搅拌桨表面标记点；由于搅拌桨是主动运动，在计算中需设置一定的转速， 因此不需要考虑流体对搅拌桨的作用力。不过当需要统计搅拌釜的操作功率时，可以通过累 加搅拌桨表面的标记点作用力与速度的乘积来确定相应的信息。

情况③：更新自由运动的颗粒信息(如颗粒的位置、速度与角速度等)；颗粒的运动与其 所受作用力和质量有关，颗粒所受作用力包括流体对颗粒的作用力以及颗粒-颗粒、颗粒-搅 拌桨或颗粒-壁面之间的相互碰撞作用力。其中，流体对颗粒的作用力可以通过累加颗粒表面 标记点的作用力得到。

优选地，搅拌反应釜内固体悬浮颗粒的模拟包括以下步骤：

流体部分在GPU端完成计算，颗粒部分在CPU端完成计算；

计算完成后，CPU端的颗粒信息传递到GPU端，进行流固耦合；

用流固耦合计算得到的欧拉网格体积力更新流体速度，然后计算得到每个颗粒的受理并 传输到CPU端以更新颗粒信息。

上述搅拌反应釜内固体悬浮颗粒的模拟步骤为一个时间步内的耦合计算过程。其中，流 固耦合主要是处理流体对固体的作用以及固体对流体的作用。

该过程用到的函数及各函数的具体功能如下：

流体部分(FluidOneStep)：计算不含体积力源项的N-S方程得到中间速度，此时的流场 就相当于不存在固体颗粒时的流场。

颗粒部分(ParticleOneStep)：该函数的功能相当于不考虑流体对颗粒作用力时单纯DEM 的计算流程。主要包括的计算过程有：邻居颗粒的搜索、计算区域的分解、颗粒网格的更新、 颗粒的受力计算、颗粒的位置以及速度等信息的更新。

流固耦合(FluidSolidInteration)：该函数需要中间速度场的信息以及颗粒的信息。由于 颗粒自由运动，我们采用直接力方法进行处理。在计算前，首先准备好颗粒表面标记点的文 件，即颗粒表面拉格朗日标记点相对于颗粒质心位置的坐标。假如计算中有不同粒径的颗粒， 那么颗粒表面标记点数目不同，需要针对不同粒径的颗粒分别生成标记点初始文件。在流固 耦合处理时，由颗粒p中心位置及第l个标记点的相对位置就可以计算得到第l个标记 点的真实位置当已知所有颗粒表面标记点的位置与预期速度，就可采用直接力 方法计算得到每个拉格朗日标记点处的体积力与周围欧拉网格的体积力源项。流体对颗粒的 作用力即为每个颗粒表面上所有标记点处的体积力之和。

更新流体速度(FluidUpdate)：通过添加体积力源项的影响，更新由函数FluidOneStep 计算得到的中间速度。

更新颗粒信息(ParticleUpdate)：在FluidSolidInteration中可以计算得到流体对颗粒的 作用力，这部分数据由GPU端传递到CPU端，并重新依据牛顿定律更新颗粒的信息。

各函数间需要传递的数据如下：

流场信息(FluidInformation)：流体三个方向的速度分量，以及各欧拉网格点的坐标；

欧拉网格体积力(EulerianBodyForce)：由标记点散布到欧拉网格上的作用力，即颗粒 对流体的作用；

颗粒信息(ParticleInformation)：需要从CPU端拷贝到GPU端，包括颗粒的位置、速度、 角速度、半径以及质量。在程序中，所有颗粒信息均放到一个结构体中。

流体对颗粒作用力(FluidForceOnParticles)：需要从GPU端拷贝到CPU端的数据，即 流体对颗粒的作用力，包括法向作用力、切向作用力与润滑力。

在浸入边界法与离散单元法耦合实现搅拌反应釜内固体悬浮颗粒的模拟中，IBM与DEM 之间需要交换的数据量较少，只需要传递颗粒的信息即可。流体计算部分由于网格数目多， 计算量较大，放在GPU上进行计算；颗粒部分由于采用界面解析的方法，颗粒数目不可能太 多(最高达Ο(10⁴))，这部分计算放在CPU端进行。因此CPU与GPU计算重叠，节省了部 分时间。

以上所述的基于浸入边界法的搅拌反应釜模拟方法，采用欧拉网格求解流动控制方程；

采用浸入边界法(ImmersedBoundaryMethod，IBM)设置边界条件；

采用大涡模拟实现搅拌反应釜内湍流流动模拟；

采用CPU-GPU异构并行计算实现搅拌反应釜的模拟；

采用浸入边界法与离散单元法(DiscreteElementMethod，DEM)的方式实现搅拌反应釜 内固体悬浮颗粒的模拟；

所述浸入边界法包括基于拉格朗日标记点的直接力方法和基于网格固含率的体积分数方 法，采用基于拉格朗日标记点的直接力方法处理运动边界(如搅拌桨表面以及悬浮颗粒表面 等)，采用基于网格固含率的体积分数方法处理静止边界(如反应釜壁面等)。

其中，边界条件为复杂边界处的边界条件，在搅拌反应釜中同时存在两种不同运动状态 的壁面，包括旋转的搅拌桨壁面和静止的搅拌釜壁面。针对不同的边界条件采用不同的浸入 边界法可以提高计算精度与效率。

采用CPU-GPU异构并行计算实现搅拌反应釜的模拟，尤其是实现搅拌反应釜的大规模 模拟。

本发明采用的浸入边界法(IBM)是在欧拉网格下求解复杂的流固耦合问题，固体边界 对流体的作用通过在控制方程中添加源项的方式实现。

控制方程为：

$\frac{\partial u}{\partial t}+u·▿u=-▿p+\frac{1}{\mathrm{Re}}{▿}^{2}u+f$

▽·u＝0

其中，f为体积力源项。

优选地，所述欧拉网格为欧拉正交网格。

优选地，所述欧拉正交网格为等距欧拉正交网格。

所述欧拉网格可采用简单的小程序生成，且所述欧拉网格在三个方向(x，y，z)的网格 宽度相等。

优选地，所述欧拉网格包含网格固含率信息。所述网格固含率中，搅拌釜器壁及其外部 区域所占据的部分为1，内部部分为0。

优选地，所述欧拉网格的网络固含率信息采用射线法进行统计。

优选地，所述欧拉网格采用交错网格时，分别统计相应变量的网格固含率，即0～1。

优选地，基于拉格朗日标记点的直接力方法中，边界条件直接施加于离散的拉格朗日标 记点处。

在基于拉格朗日标记点的直接力方法中，如果边界是运动的，那么每一时间步均需要更 新标记点的信息。因为只有固体的表面分布标记点，所以相对于欧拉网格来说，拉格朗日标 记点的数目较小。这种方法能够十分有效的抑制运动边界计算中出现的非物理的力的震荡， 但是对于静止边界则没有类似的力的震荡。

另外，当浸入边界十分靠近整个计算区域的边界时，拉格朗日标记点与计算区域边界的 距离小于狄拉克函数的支持宽度时，就会出现没有足够的欧拉网格用于速度的差值以及体积 力的散布。这就需要作特殊的处理，但势必会引入不必要的误差。因而需要采取基于网格固 含率的体积分数方法。

优选地，基于网格固含率的体积分数方法中体积力的计算通过浸入边界处固体的预期速 度计算得到。

优选地，基于网格固含率的体积分数方法中流体与固体区域的速度通过网格固含率来进 行平滑过渡。

当固体区域静止不动时，基于网格固含率的体积分数方法中每个网格的固含率在计算中 也保持不变，因而网格固含率也不需要每个时间步均重新计算。

基于拉格朗日标记点的直接力方法与基于网格固含率的体积分数方法分别更加适合于处 理运动边界与静止边界。在混合浸入边界法中，欧拉网格上的体积力等于运动边界的作用力 分量加上静止边界的作用力分量。

优选地，基于拉格朗日标记点的直接力方法中拉格朗日标记点的生成方法包括以下步骤：

(1)生成搅拌桨的三维构体并导出相应的结构信息；其中，网格信息是指用软件生成的 存储有三维构体信息的文件。

(2)利用导出的结构信息生成搅拌桨表面网格，并导出搅拌桨表面网格的位置坐标，然 后计算得到拉格朗日标记点的控制体积。

在浸入边界法中，尽量均匀的标记点能够得到更为精确的计算结果，但是对于复杂的几 何表面，若要生成均匀的标记点几乎是不可能的，因此对于几何形状种类很多的搅拌桨来说， 可以借助软件来生成标记点，如SolidWorks、Pro/E与Gambit等。首先由SolidWorks或Pro/E 生成三维构体，导出STL格式文件，再将生成STL文件导入到Gambit软件中，进而生成搅 拌桨表面的网格。生成的网格的每个格点间的距离大致等于欧拉网格的宽度，导出搅拌桨表 面网格的位置坐标即为拉格朗日标记点的初始坐标，拉格朗日标记点的控制体积由以下公式 计算得到：

控制体积＝搅拌桨表面积/标记点的数目×欧拉网格宽度

优选地，所述大涡模拟中采用拉格朗日动态亚格子模型。

优选地，所述大涡模拟中需添加大涡模拟-壁函数。

通常情况下，搅拌反应釜的工作雷诺数(Re)较高，反应釜内流动尤其是搅拌桨区域附 近均达到湍流状态。现今，最精确的模拟湍流流动的方法是直接数值模拟(DirectNumerical Simulation，DNS)，它能够分辨湍流各个尺度的特征，但其总的计算量与Re³成正比，对于 当前的计算机条件还很难实现较高Re的湍流直接数值模拟。而雷诺平均模型(Reynolds AverageNavier-Stokes，RANS)虽然计算量较小，但是计算精度太差，无法精确捕捉搅拌反 应釜内湍流的细节。

本发明采用的大涡模拟(LargeEddySimulation，LES)计算量介于DNS与RANS之间， 其基本思想是直接计算大尺度脉动，而小尺度脉动的影响通过模型来考虑。LES被证明是一 种处理搅拌反应釜湍流流动模拟的有效方法。要实现LES需要构建亚格子应力模型计算得到 湍流粘度，最常用的亚格子应力模型是Smagorinsky提出的涡粘模型，但其无法反应湍流特 征在近壁处的渐进特性，会造成耗散过大。

本发明在IBM框架下实现LES，由于搅拌釜内流动问题比较复杂，难以找到统计均匀的 空间方向，因而采用拉格朗日动态亚格子模型，该模型可以沿质点轨迹平均进而动态适应局 部湍流结构并由此确定模型系数。

另外为了减小壁面附近的网格分辨率，还需要添加大涡模拟-壁函数。

对于边界层类型的流动，粘性壁面区域非常重要，其会对外流区产生极大的影响。如果 解析近壁区域，大涡模拟的网格数将与直接数值模拟处于同一量级，这就失去了大涡模拟的 优越性。为了保持大涡模拟的优势，并尽量减少计算量，在近壁区域常采用近似边界条件模 化壁面层对外流的影响。本发明采用幂函数来模化近壁速度。为了添加壁函数，壁面处不是 采用无滑移边界，而是采用滑移边界条件，这是为了在浸入边界处给定合适的速度使其能够 为外部大涡模拟的计算提供合适的壁面剪切应力。在浸入边界法中施加滑移边界只需修改拉 格朗日标记点处的预期速度。

优选地，所述CPU-GPU异构并行计算包括：

(1)流体求解的GPU化；

(2)流固耦合的GPU化；

(3)流体求解的并行处理；

(4)流固耦合的并行处理；

其中，流体求解的GPU化主要是采用CUDA在GPU上进行计算；为了减少GPU与CPU 之间的数据拷贝，流固耦合部分也在GPU上计算。

复杂流动横跨不同的空间与时间尺度，即便采用LES计算量依然十分巨大。另外，若要 实现颗粒流体两相流的界面解析模拟，单独采用串行的CPU程序计算规模有限。故本发明采 用CPU-GPU异构并行方程来实现搅拌反应釜内复杂流动的大规模模拟。

搅拌反应釜内，整个流固耦合过程包含流体求解部分、固体部分以及流固耦合部分，由 于几部分间算法差异明显，其GPU化以及并行算法需要分别进行设计。为了减少数据在CPU 与GPU之间拷贝的消耗，流体的求解计算均在GPU上进行。

优选地，流体求解的GPU化中GPU的每一个线程负责处理一个欧拉正交网格单元。

其中，不同的变量均按照相同的顺序连续存储。GPU程序需要被分割成不同的程序段， 以便在CPU端实现全局同步。且在GPU端的程序运行过程中，CPU主要是担任协调控制的 作用，同时也会执行与GPU核并行的计算任务。

为了实现大规模数值模拟，需要借助于高性能并行计算技术，同样的对于流体求解与流 固耦合采用不同的并行处理方式。

优选地，流体求解的并行处理中采用规则的矩形计算区域。

优选地，流体求解的并行处理中将整个计算区域划分成多个亚计算区域，每个亚计算区 域为一个进程；其中“多个”为2个或2个以上，其可以由本领域技术人员根据经验和实际 情况确定，在本领域是清楚的。

优选地，流体求解的并行处理中每个进程运行在一个CPU以及相应的GPU上。

优选地，流体求解的并行处理中相邻进程间通过MPI(MultiPointInterface，多点接口) 进行数据传输。

优选地，流固耦合的并行处理中采用了两种浸入边界法。

优选地，两种浸入边界法的并行方式不同。

优选地，流固耦合的并行处理中采用的两种浸入边界法为基于拉格朗日标记点的直接力 方法和基于网格固含率的体积分数方法。

优选地，基于网格固含率的体积分数方法中网格固含率与欧拉网格存储在一起。

所述网格固含率与欧拉网格存储在一起，即网格固含率存在于欧拉网格的每个网格位置。 并且每个网格的网格固含率与周围的网格间没有关联，因此体积分数方法的并行处理没有任 何额外的难度，其GPU计算与并行处理与流体计算一致。

优选地，基于拉格朗日标记点的直接力方法中采用部分速度插值和体积力散布的方式来 进行浸入边界法的并行处理。

所述部分速度插值和体积力散布均以标记点进行并行。

在体积力散布过程中有可能会出现不同的线程向全局存储器的同一位置进行写操作，这 种情况可以通过对全局变量的原子操作实现对共享可写数据的互斥操作。

基于拉格朗日标记点的直接力方法由于存在自由运动的拉格朗日标记点，因此增加了并 行处理的复杂度；另外，欧拉与拉格朗日变量之间需要进行交换，每个拉格朗日标记点与周 围多个欧拉网格相互关联，同时每个欧拉网格也与多个拉格朗日标记点存在关系。因而需要 采用部分速度插值和体积力散布的方式来进行简化。该方法能够较为简便的处理浸入边界穿 越进程边界的问题而不增加通信量。

搅拌反应釜内的颗粒悬浮在化工过程中有广泛的应用。液固两相流的颗粒解析模拟不仅 能够记录每个颗粒的具体信息，还能够捕捉到颗粒周围流场的细节。因为需要求解颗粒间隙 的流场，欧拉网格的宽度一般比颗粒直径小一个数量级。

本发明中采用IBM-DEM耦合的方式处理搅拌釜固体悬浮的颗粒解析模拟。

优选地，采用浸入边界法与离散单元法耦合的方式实现搅拌反应釜内固体悬浮颗粒的模 拟中，浸入边界法直接在自由运动颗粒表面施加边界条件。

优选地，采用浸入边界法与离散单元法的方式实现搅拌反应釜内固体悬浮颗粒的模拟中， 离散单元法处理颗粒的运动与碰撞。

与现有技术相比，本发明具有以下有益效果：

(1)本发明所述方法的网格生成简便并且改变桨型不需要重新生成网格。与传统的MRF 或者SM方法不同的是，其只需要重新得到搅拌桨拉格朗日标记点的坐标即可，使工作量大 大降低。

(2)本发明所述方法具有可以处理任意复杂运动边界的能力。本发明采用浸入边界方法 可以直接在复杂边界处施加无滑移边界条件，并且搅拌桨几何结构的复杂程序并不会增加实 现的困难，可以完全避免MRF与SM在精确模拟强烈的转子定子相互作用或者是当搅拌桨与 其他静止附件有相互嵌入时所存在的缺陷。

(3)本发明所述混合浸入边界方法，即采用基于网格固含率的体积分数方法处理静止边 界，采用基于拉格朗日标记点的直接力方法处理运动边界，不仅提高了计算效率而且保证了 计算精度。在直接力方法中，边界条件直接在离散的拉格朗日标记点处施加。如果边界是运 动的，那么每一时间步均需要更新标记点的信息。因为只有固体的表面分布标记点，所以相 对于欧拉网格来说，标记点数目较小。这种方法能十分有效的抑制动边界计算中出现的非物 理的力的震荡。在体积分数方法中，体积力的计算通过浸入边界处固体的预期速度直接计算 得到，流体与固体区域的速度通过网格固含率来实现平滑过渡。当固体区域静止不动时，每 个网格的固含率在计算中也保持不变，因而网格固含率也不需要每个时间步均重新计算。

(4)本发明采用CPU-GPU大规模并行计算，使得工业尺度搅拌釜的模拟成为可能。传 统情况下，常采用CPU并行计算且由于贴体网格每个网格点的计算操作要远多于欧拉正交网 格，所以计算网格数不会太多，一般为10⁶量级。而本发明采用CPU-GPU并行计算后可以计 算的网格数可达到10⁸～10⁹量级，因此可以实现工业尺度搅拌反应釜的模拟。

(5)本发明搅拌釜内颗粒悬浮的界面解析模拟中，颗粒表面的边界条件直接采用浸入边 界法施加。相对于传统的液固两相流的欧拉-欧拉模型以及欧拉-拉格朗日模型，其不仅能够 记录每个颗粒的具体信息，并且能够捕捉到每个颗粒周围流场的细节。流固之间的相互作用 通过在固体颗粒表面直接施加边界条件，不需要借助于经验模型。

(6)本发明提供的基于浸入边界法的搅拌反应釜的模拟方法，结合高性能计算能够实现 工业尺度的搅拌反应釜的模拟；搅拌釜数值模拟能够得到反应釜内部详细的流场与温度场信 息，能够为反应釜的优化设计提供依据。

附图说明

图1是采用传统多重网格法与滑移网格法进行搅拌反应釜模拟的计算区域示意图；

图2是两种不同浸入边界法设置边界条件示意图，其中(a)为基于拉格朗日标记点的直 接力方法，(b)为基于网格固含率的体积分数方法；

图3是本发明所述方法中准备阶段和数值计算阶段计算流程图；

图4是浸入边界法进行搅拌桨与反应釜模拟的计算结果图，其中case1：搅拌桨与反应釜 壁均采用体积分数方法，case2：搅拌桨与反应釜壁均采用直接力方法，case3：搅拌桨与反 应釜壁分别采用直接力方法与体积分数方法；

图5是GPU并行计算示意图；

图6是CPU与GPU性能对比图；

图7是二维笛卡尔网格并行计算区域划分示意图；

图8是采用浸入边界法与离散单元法耦合的方式实现搅拌反应釜内固体悬浮颗粒的模拟 示意图；

图9是采用大涡模拟计算的Re＝7300的Rushton桨搅拌槽中搅拌桨横截面处的瞬时涡量 图；

图10是采用大涡模拟计算的Re＝7300的Rushton桨搅拌槽中通过转轴横截面处的瞬时涡 量图；

图11是采用不同方法模拟得到的Re＝7300的Rushton搅拌槽的相平均速度，其中(a)、 (b)和(c)为径向速度分量在三个不同径向位置处(r/T＝0.183，r/T＝0.25，r/T＝0.317)沿轴 向的分布，(d)、(e)和(f)为切向速度分量在三个不同径向位置处(r/T＝0.183，r/T＝0.25， r/T＝0.317)沿轴向的分布；

图12是工业尺度搅拌釜内瞬时涡量图。

具体实施方式

下面结合附图并通过具体实施方式来进一步说明本发明的技术方案。本领域技术人员应 该明了，所述实施例仅仅用于帮助理解本发明，不应视为对本发明的具体限制。

传统多重网格法与滑移网格法进行搅拌反应釜模拟的计算区域示意图，在这两种方法 中，计算区域被划分成静止区域与运动区域，计算结果在内外部界面处直接进入匹配，严重 影响了计算效率以及计算精度。另外如果需要精确的模拟强烈的转子定子相互作用或者是当 搅拌桨与其他静止附件有相互嵌入时，这些方法不适用。

本发明采用的浸入边界法(IBM)在欧拉网格下求解复杂的流固耦合问题，固体边界对 流体的作用通过在控制方程中添加源项的方式实现。

控制方程为：

$\frac{\partial u}{\partial t}+u·▿u=-▿p+\frac{1}{\mathrm{Re}}{▿}^{2}u+f$

▽·u＝0

其中，f为体积力源项。

以下通过实施例来说明本发明所述的基于浸入边界法的搅拌反应釜模拟方法。

实施例1：

1、准备阶段：

计算参数的测定和设置，计算参数包括数值参数(时间步长、网格数以及进程数等)与 物理参数(搅拌釜结构参数、流体的密度、流体的粘度、颗粒直径和颗粒密度等)。

准确的搅拌釜结构结构参数可以从其设计图纸中获得；在没有图纸的情况下，我们可以 采用长度测量工具，如千分尺等进行长度的测量。主要的结构参数有搅拌反应釜的高度与直 径，搅拌桨的几何形状与尺寸等。

流体的粘度采用毛细管粘度计、落球黏度计或同轴圆筒粘度计中任意一种进行测量。

流体的密度采用密度计进行测量。

颗粒直径采用激光粒度仪进行测量。

颗粒密度采用粉末和微小颗粒测试领域的密度仪进行测量。

程序计算区域为矩形区域，六个面可设置为无滑移边界，有时顶部平面应设置滑移壁 面，具体视情况而定。

生成流体网格文件，其中包含每个网格点的位置坐标，速度与压力初始数值，网格固含 率数据。

由于网格为等距欧拉正交网格，可以采用很简单的小程序生成，其中网格固含率采用射 线法进行统计。当采用并行计算时，需要把计算网格划分成较小的亚计算区域，具体划分视 情况而定。

生成拉格朗日标记点信息，其包括搅拌桨表面的标记点信息，另外在搅拌釜内颗粒悬浮 模拟时，需给出颗粒表面标记点的信息；拉格朗日标点信息生成的具体步骤如下：

(1)借助SolidWorks或Pro/E等商业软件生成搅拌桨的三维构体图，并导出STL格式； 对于搅拌反应釜的静止部件(如器壁与挡板等)以及搅拌反应釜的运动部件(如搅拌桨等) 需要分别构图，然后导出相应的STL文件。

(2)对于静止部件，采用单独的程序来处理STL文件，假设在实际计算中其中一个方 向的网格数为N，则在处理STL文件时，我们在此方向上等距离离散成m×N个格点并采用 射线法来确定格点是位于固体区域还是流体区域，其中位于固体区域内的格点标记为1，而 位于流体区域内的格点标记为0。那么可以统计位于一个欧拉网格内标记为1的格点数目为 n，因此可以得到该欧拉网格的固含率为n/m³。

对于运动部件，将生成的STL导入Gambit中，生成三角形网格。网格结点的坐标即为 拉格朗日标记点的坐标，需尽量保证拉格朗日点之间的间距大致与欧拉网格宽度相等。

计算得到的欧拉网格固含率存储在初始欧拉网格中，搅拌桨表面拉格朗日标记点存储在 另外一个单独的文件中，程序在计算开始时读取这些初始数据用于计算。

准备计算资源，包括CPU资源和GPU资源，在并行计算中要求有足够的CPU以及GPU 计算机资源。

2、数值计算：

(1)读取准备生成的流体网格文件和计算参数，并将其由CPU端传递到GPU端；

(2)在不考虑控制方程中体积力源项的情况下求解纳维-斯托克斯方程得到预估速度。

本发明中N-S方程的求解采用分数步方法，N-S方程的粘性项采用四阶中心差分格式， 计算区域边界处降为两阶；N-S方程的对流项采用三阶迎风格式，在边界处采用二阶中心差 分格式；控制方程的时间推进采用显示三阶Runge-Kutta方法。

在并行计算中，亚计算区域中需要设置两层虚拟网格用于存储邻居进程的数据，该部分 数据需要利用MPI进行传输。

(3)求解压力泊松方程得到压力修正值，进一步更新预估速度；

求解压力泊松方程采用二阶中心差分格式离散，并行计算中每一个迭代步内均需要更新 虚拟网格数据。

(4)流固耦合处理，即考虑搅拌桨、反应釜壁面以及内部颗粒等对留流场的影响。采用 直接力方法处理运动边界得到作用力f₁，采用体积分数方法处理静止边界得到流固界面处的 体积力f₂，将f₁与f₂累加得到作用力f并更新流场速度；

两种不同浸入边界法设置边界条件如图2所示。

采用直接力方法处理运动边界包括以下步骤：

首先，通过狄拉克函数把欧拉网格上的速度u转换到拉格朗日标记点处的U；然后计算 拉格朗日标记点处U的作用力；再将拉格朗日标记处的作用力散布到欧拉网格上得到作用力 f₁；

其中，拉格朗日标记点处的U为拉格朗日标记点处U 的作用力是通过预期速度V^d与当前计算速度的差值计算得到，即其 中RHS包含了对流项、压力项以及粘性项；将拉格朗日标记处的作用力散布到欧拉正交网 格上得到作用力f₁，$f_1\left(x\right)=\underset{l=1}{\overset{N_l}{\Sigma }}F\left(X_l\right){\delta }_h(x-X_l){\mathrm{\Delta V}}_l.$

采用体积分数方法处理静止边界是通过网格固含率插值得到流固界面处的体积力f₂；其 中，f₂＝α(u_d-u)/Δt。

(5)通过作用力f和流场速度更新浸入边界的位置与速度；包括以下三种情况中任意一 种或至少两种：

情况①：搅拌反应釜避免始终保持静止时，不更新网格固含率，此部分不需要作任何处 理；

情况②：更新搅拌桨表面标记点；由于搅拌桨是主动运动，在计算中需设置一定的转 速，因此不需要考虑流体对搅拌桨的作用力。不过当需要统计搅拌釜的操作功率时，可以通 过累加搅拌桨表面的标记点作用力与速度的乘积来确定相应的信息。

情况③：更新自由运动的颗粒信息；颗粒的运动与其所受作用力和质量有关，颗粒所受 作用力包括流体对颗粒的作用力以及颗粒-颗粒、颗粒-搅拌桨或颗粒-壁面之间的相互碰撞作 用力。其中，流体对颗粒的作用力可以通过累加颗粒表面标记点的作用力得到。

以上准备阶段和数值处理阶段的流程图如图3所示。

本发明采用两种不同浸入边界法处理边界问题的计算结果如图4所示，同时，本发明还 对搅拌桨与反应釜壁均采用体积分数方法和搅拌桨与反应釜壁均采用直接力方法两种情况进 行了计算，结果如图4所示。

可以看出采用不同的浸入边界法分别处理体系中的运动边界与静止边界，可以提高计算 效率，同时也能够保证计算精度，其结果在三个算例中是最好的。

同时，采用不同的浸入边界法分别处理体系中的运动边界与静止边界其在搅拌釜模拟中 的效率优于另外两种方法，其结果如表1所示。表格中的数据是每个算例运行100个时间步 所需时间，从表1中可以看出采用混合浸入边界法的计算时间大大降低。

表1：不同浸入边界法在搅拌釜模拟中的效率对比表

3、将数值计算阶段得到浸入边界的位置与速度应用于搅拌反应釜内湍流流动模拟、搅 拌反应釜的大规模模拟或搅拌反应釜内固体悬浮颗粒的模拟中任意一种模拟或至少两种的模 拟中。

(1)将数值计算阶段得到浸入边界的位置与速度应用于搅拌反应釜内湍流流动模拟中：

本发明中采用拉格朗日动态亚格子模型，该模型是一种沿质点轨迹作平均从而动态确定 模型系数的方法。该方法不需要自定义的参数，亚格子模型系数可以由可解速度场计算得 到；其在近壁区内或者层流区内可自动调整或者关闭，不需要衰减函数。在近壁区域中，本 发明采用幂函数来模化近壁速度。为了添加壁函数，壁面处不是采用无滑移边界，而是采用 滑移边界条件，这是为了在浸入边界处给定合适的速度使得能够为外部大涡模拟的计算提供 合适的壁面剪切应力。在浸入边界法中施加滑移边界可以较为简单的通过修改拉格朗日标记 点处的预期速度即可。

本发明采用大涡模拟对Re＝7300的Rushton桨搅拌反应釜进行湍流流动模拟，其结果如 图9、图10、图11和图12所示，从图中可以看出采用大涡模拟能够得到搅拌桨附近详细的 湍流流动，能够捕捉到不同尺度的涡结构。另外从图12中计算结果的定量分析可以看出， 采用大涡模拟的平均速度与实验结果更加吻合，说明了大涡模拟在湍流流动模拟中的能力。

(2)将数值计算阶段得到浸入边界的位置与速度应用于搅拌反应釜的大规模模拟，本发 明采用CPU-GPU异构并行计算实现搅拌反应釜的模拟，其并行计算如图5所示。CPU-GPU 异构并行计算相对于常规的CPU和GPU，可以实现工业尺度搅拌反应釜的模拟，CPU和 GPU计算对比如图6所示，从图中可以看出采用GPU后计算时间大大降低，GPU相对于单 核CPU的加速比为7～20倍，相对于CPU四核的加速比也能够达到5倍左右。网格数量越 多，加速效果越明显。

并行计算-流体求解并行处理：

此过程中流体求解采用规则的矩形计算区域，并被均匀划分到亚计算区域中，如图7所 示。每一个亚计算区域为一个进程，并运行在一个CPU以及相应的GPU上。每个亚计算区 域边界处的网格计算均需要邻近的网格作差分运算，网格层数由采用的数值差分方法决定， 所采用二阶中心差分格式需要一层虚拟网格，而采用四阶中心差分则需要两层虚拟网格。在 每一次迭代中，虚拟网格中的数据均需要从相邻进程中通过MPI传输。由于流体计算区域均 匀划分，拉格朗日标记点在各进程间可能分布不均，但由于流体的求解占据了绝大部分的时 间，所以不会造成严重的负载不均衡的问题。

并行计算-浸入边界法并行处理：

本发明采用基于网格固含率的体积分数方法处理静止边界，采用基于拉格朗日标记点的 直接力方法处理运动边界。在体积分数方法中，网格固含率存在于每个网格位置，并与周围 的网格间没有关联，所以这种方法的并行处理没有任何额外的难度。对于基于拉格朗日标记 点的直接力方法中，采用一种不完全速度插值与分布力的方式。在这种实现方式中只需要一 层虚拟网格，速度的插值与力的散布均需作相应的修改。当前采用的不完全速度插值与体积 力散布的方式在并行过程中不会引入额外的通信，从而使得并行更加高效并且更易于实现。

(3)将数值计算阶段得到浸入边界的位置与速度应用于搅拌反应釜内固体悬浮颗粒的模 拟中。

本发明采用浸入边界法与离散单元法耦合的方式实现搅拌反应釜内固体悬浮颗粒的模 拟，模拟流程如图8所示。浸入边界法与离散单元法耦合两部分之间需要交换的数据量较 少，只需要传递颗粒的信息即可。流体计算部分由于网格数目多，计算量较大，放在GPU 上进行计算；颗粒部分由于采用界面解析的方法，颗粒数目不可能太多(最高达O(10⁴))， 这部分计算放在CPU端计算。当流体部分(FluidOneStep)与颗粒部分(ParticleOnestep) 均计算完成后，CPU端的颗粒信息，包括颗粒的位置、速度、角速度、半径与质量等信息传 递到GPU端，然后实现流固耦合部分(FluidSolidInteration)。流固耦合主要是处理流体对 固体的作用以及固体对流体的作用，其中计算得到的欧拉网格体积力用来更新流体速度 (FluidUpdate)，计算得到的每个颗粒的受力需要传输到CPU端用于更新颗粒信息(Particle Update)。至此就完成了一个时间步内的耦合计算过程。

流体部分(FluidOneStep)：计算不含体积力源项的N-S方程得到中间速度，此时的流 场就相当于不存在固体颗粒时的流场。

颗粒部分(ParticleOneStep)：该函数的功能相当于不考虑流体对颗粒作用力时单纯 DEM的计算流程。主要包括的计算过程有：邻居颗粒的搜索、计算区域的分解、颗粒网格 的更新、颗粒的受力计算、颗粒的位置以及速度等信息的更新。

流固耦合(FluidSolidInteration)：该函数需要中间速度场的信息以及颗粒的信息。由于 颗粒自由运动，我们采用直接力方法进行处理。在计算前，首先准备好颗粒表面标记点的文 件，即颗粒表面拉格朗日标记点相对于颗粒质心位置的坐标。假如计算中有不同粒径的颗 粒，那么颗粒表面标记点数目不同，需要针对不同粒径的颗粒分别生成标记点初始文件。在 流固耦合处理时，由颗粒p中心位置及第l个标记点的相对位置就可以计算得到第l个 标记点的真实位置当已知所有颗粒表面标记点的位置与预期速度，就可采用直 接力方法计算得到每个拉格朗日标记点处的体积力与周围欧拉网格的体积力源项。流体对颗 粒的作用力即为每个颗粒表面上所有标记点处的体积力之和。

更新流体速度(FluidUpdate)：通过添加体积力源项的影响，更新由函数Fluid_One_Step 计算得到的中间速度。

更新颗粒信息(ParticleUpdate)：在FluidSolidInteration中可以计算得到流体对颗粒的 作用力，这部分数据由GPU端传递到CPU端，并重新依据牛顿定律更新颗粒的信息。

各函数间需要传递的数据如下：

流场信息(FluidInformation)：流体三个方向的速度分量，以及各欧拉网格点的坐标；

欧拉网格体积力(EulerianBodyForce)：由标记点散布到欧拉网格上的作用力，即颗粒 对流体的作用；

颗粒信息(ParticleInformation)：需要从CPU端拷贝到GPU端，包括颗粒的位置、速 度、角速度、半径以及质量。在程序中，所有颗粒信息均放到一个结构体中。

流体对颗粒作用力(FluidForceOnParticles)：需要从GPU端拷贝到CPU端的数据， 即流体对颗粒的作用力，包括法向作用力、切向作用力与润滑力。

(4)后处理阶段，模拟完成后输出每个欧拉网格点的信息以及每个颗粒的信息。

后处理可以采用Tecplot、Paraview以及Origin等工具进行处理，同时也可以写一些小的 程序用于辅助数据的处理。

本发明模拟的工业尺度搅拌釜内瞬时涡量如图12所示，左侧两个插图同时显示了垂直 平面内不同位置处的流线图，从图中可以看出清晰的看出搅拌桨附近区域的涡量明显大于其 他区域。在每一个搅拌桨后部均形成了两个较大的尾涡。由于换热蛇管的阻挡作用，蛇管与 搅拌桨间流动区域的流速远远大于蛇管与侧壁面之间区域的流速。我们注意到中层搅拌桨排 出流上下方的速度值很小，这可能是由于流体粘度较大或者计算时间过短导致的。换热蛇管 后面的狭小区域内的流动细节，在DNS中仍然能够很好的捕捉到。从图中的插图可以看 出，由于流体流经盘管形成绕流，在盘管后部形成明显的尾涡。从该算例可以看出，采用混 合浸入边界方法能够很容易处理含有复杂运动与静止边界的流动与传热现象。

综上所述，可以看出本发明所述方法的网格生成简便并且改变桨型不需要重新生成网 格，其只需要重新得到搅拌桨拉格朗日标记点的坐标即可，使工作量大大降低；本发明采用 浸入边界方法可以直接在复杂边界处施加无滑移边界条件，并且搅拌桨几何结构的复杂程序 并不会增加实现的困难，可以完全避免MRF与SM在精确模拟强烈的转子定子相互作用或者 是当搅拌桨与其他静止附件有相互嵌入时所存在的缺陷；本发明所述混合浸入边界方法，即 采用基于网格固含率的体积分数方法处理静止边界，采用基于拉格朗日标记点的直接力方法 处理运动边界，不仅提高了计算效率而且保证了计算精度；本发明采用CPU-GPU大规模并 行计算，使得工业尺度搅拌釜的模拟成为可能；本发明搅拌釜内颗粒悬浮的界面解析模拟 中，不仅能够记录每个颗粒的具体信息，并且能够捕捉到每个颗粒周围流场的细节。流固之 间的相互作用通过在固体颗粒表面直接施加边界条件，不需要借助于经验模型。

申请人声明，本发明通过上述实施例来说明本发明的详细工艺设备和工艺流程，但本发 明并不局限于上述详细工艺设备和工艺流程，即不意味着本发明必须依赖上述详细工艺设备 和工艺流程才能实施。所属技术领域的技术人员应该明了，对本发明的任何改进，对本发明 产品各原料的等效替换及辅助成分的添加、具体方式的选择等，均落在本发明的保护范围和 公开范围之内。