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

Abstract

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

Description

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

技术领域

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

背景技术

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

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

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

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

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

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

Verzicco等人(FlowinanImpeller-StirredTankUsinganImmersed-BoundaryMethod.AicheJ，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-SlipFlowBoundarywithanExternalForce-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}{\Σ}}F\left(X_l\right){\mathrm{\δ}}_h(x-X_l){\mathrm{\Δ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{\∂u}{\∂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)较高，反应釜内流动尤其是搅拌桨区域附近均达到湍流状态。现今，最精确的模拟湍流流动的方法是直接数值模拟(DirectNumericalSimulation，DNS)，它能够分辨湍流各个尺度的特征，但其总的计算量与Re³成正比，对于当前的计算机条件还很难实现较高Re的湍流直接数值模拟。而雷诺平均模型(ReynoldsAverageNavier-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{\∂u}{\∂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}{\Σ}}F\left(X_l\right){\mathrm{\δ}}_h(x-X_l){\mathrm{\Δ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端用于更新颗粒信息(ParticleUpdate)。至此就完成了一个时间步内的耦合计算过程。

流体部分(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大规模并行计算，使得工业尺度搅拌釜的模拟成为可能；本发明搅拌釜内颗粒悬浮的界面解析模拟中，不仅能够记录每个颗粒的具体信息，并且能够捕捉到每个颗粒周围流场的细节。流固之间的相互作用通过在固体颗粒表面直接施加边界条件，不需要借助于经验模型。

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

Claims (10)

1.一种基于浸入边界法的搅拌反应釜模拟方法，其特征在于，所述方法包括以下步骤：

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

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

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

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

2.根据权利要求1所述的搅拌反应釜模拟方法，其特征在于，步骤(1)中计算参数包括数值参数和物理参数；

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

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

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

优选地，搅拌釜结构参数采用长度测量工具测定；

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

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

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

优选地，颗粒密度采用密度仪进行测量；

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

优选地，所述流体网格为欧拉网格；

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

优选地，步骤(1)中拉格朗日标记点信息还包括颗粒表面标记点的信息；

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

优选地，步骤(1)中准备计算资源包括CPU资源和GPU资源。

3.根据权利要求1或2所述的搅拌反应釜模拟方法，其特征在于，步骤(2)数值计算包括以下步骤：

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

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

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

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

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

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

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

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

4.根据权利要求1-3任一项所述的搅拌反应釜模拟方法，其特征在于，搅拌反应釜内固体悬浮颗粒的模拟包括以下步骤：

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

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

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

优选地，CPU端的颗粒信息包括颗粒的位置、速度、角速度、半径和质量。

5.根据权利要求1-4任一项所述的搅拌反应釜模拟方法，其特征在于，所述方法采用欧拉网格求解流动控制方程；

采用浸入边界法设置边界条件；

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

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

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

所述浸入边界法包括基于拉格朗日标记点的直接力方法和基于网格固含率的体积分数方法，采用基于拉格朗日标记点的直接力方法处理运动边界，采用基于网格固含率的体积分数方法处理静止边界。

6.根据权利要求1-5任一项所述的搅拌反应釜模拟方法，其特征在于，所述欧拉网格为欧拉正交网格；

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

优选地，所述欧拉网格包含网格固含率信息；

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

优选地，所述欧拉网格采用交错网格时，分别统计变量(u，v，w)的网格固含率。

7.根据权利要求1-6任一项所述的搅拌反应釜模拟方法，其特征在于，所述基于拉格朗日标记点的直接力方法中，边界条件直接施加于离散的拉格朗日标记点处；

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

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

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

(1)生成搅拌桨的三维构体并导出相应的结构信息；

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

8.根据权利要求1-7任一项所述的搅拌反应釜模拟方法，其特征在于，所述大涡模拟中采用拉格朗日动态亚格子模型；

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

9.根据权利要求1-8任一项所述的搅拌反应釜模拟方法，其特征在于，所述CPU-GPU异构并行计算包括：

(1)流体求解的GPU化；

(2)流固耦合的GPU化；

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

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

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

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

优选地，流体求解的并行处理中将整个计算区域划分成多个亚计算区域，每个亚计算区域为一个进程；

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

优选地，流体求解的并行处理中相邻进程间通过MPI进行数据传输；

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

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

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

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

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

10.根据权利要求1-9任一项所述的搅拌反应釜模拟方法，其特征在于，采用浸入边界法与离散单元法耦合的方式实现搅拌反应釜内固体悬浮颗粒的模拟中，浸入边界法直接在自由运动颗粒表面施加边界条件；

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