首页/文章/ 详情

界面追踪方法Level Set与VOF在气泡流动模拟的效果比较

4月前浏览5231

对于两相流模拟,模型主要分为两大类:高相分数模型和界面捕捉类模型。当我们关注水中的含气量(气泡界面及气泡形状可忽略),则采用高相分数模型,此模型适用于气泡特别多的流动问题。对于有明确边界的流体-流体问题,基本需要考虑如何捕捉边界面。常用的界面捕捉模型包括LS(Level Set)方法和VOF(Volume of Fluid)方法。

多相流模拟软件,首先就是针对此类有边界面的问题。目前主流的商业CFD软件大多采用VOF方法,而定位于多相流仿真的国产通用流体仿真软件Virtualflow采用Level Set方法进行界面流仿真。

1. Level Set方法    

Level Set方法是基于空间曲面的隐函数表达。

在LS方法中,每一个时间步都要重新初始化LS方程,在时刻tn 求得的LS函数与控制方程一起求解得到下一时刻的LS函数,这些初始化的过程中总伴随着界面位置的移动,会造成质量损失,导致质量不守恒。而改善初始化步骤来矫正质量守恒又会增加计算时间,提升计算成本。同时,因为LS方法采用的是光滑的距离函数来捕捉相界面,各个物理量可以在界面上光滑连续地过渡,且相界面的捕捉效果好。

2. VOF方法    

在VOF方法中,用来划分两相界面的函数是体积分数α,表示的是单个网格内的液体体积与这个网格总体积的比值。

若求出整个网格相分数,如图1(a)来构造界面,会发现体积分数α在空间上是一个阶梯函数,在空间上是不连续的,从而重构出来的相界面,如图1(b)是间断的,两个相邻网格的界面是不连续的,且物理量在通过界面时也是不连续的,这个现象称为寄生流动,目前VOF方法的主要工作就是缓解数值方法造成的寄生流动现象。


图1 相分数空间分布(a)及其界面重构(b)

与LS方法类似,根据相分数可以得到界面上的单位法向量和曲率以及计算域中的密度和粘度。

为了解决物理量在界面两端不连续的问题,引入连续表面力模型,通过以下公式将界面内的压力表示为压力的连续函数。

式中c是界面处的位置函数(图2),其表示为:

图2 CSF模型下c位置函数(左)及连续压力函数(右)

由以上公式可以推出曲面微元上的表面张力。

此外,在OpenFOAM中,为了求解动量守恒方程中的压力项和体积力项,定义prgh,如下:

式中为网格位置矢量,对该公式求其梯度得到:

该公式可以直接带入动量守恒方程中进行计算。

在OpenFOAM中,使用VOF方法后在控制方程中添加了一个求解α的相方程:

为了界面的尖锐,OpenFOAM采用Waller提出的方法,在相方程中添加人工对流项,从而保证界面的清晰。

其中:

c为压缩因子,值为0时表示不存在人工压缩,给c赋值后有利于提高界面清晰度,但同时也会提高计算成本和产生收敛问题。

3. 界面捕捉效果:LS vs VOF    

本文主要讨论通用流体仿真软件Virtualflow中用到的LS方法和开源软件OpenFOAM中用到的VOF方法。为了验证LS方法和VOF方法对界面捕捉的效果,下面展示Albadawi文献中采用这两种方法计算模拟的气泡变化过程,并与实验进行了对比分析。在计算域底部中心一个小孔以恒定体积流率喷射气泡,由于压力、浮力和表面张力的共同作用,气泡会经历产生,变形和分离的过程。计算域及其物性参数如下:

图3 计算域及其边界条件

表2 两相流物性参数

取tDet 为气泡分离时间,t/tDet =0 为初始时间,此时气泡已经是轴对称的球形,实验中各个时刻气泡形状如图4所示。

图4  不同时刻气泡形状

图5  LS方法(TransAT软件)模拟结果与实验数据对比

图6   VOF方法(OpenFOAM)模拟结果与实验数据对比

由图5与图6可以看出:与实验结果对比,LS方法对界面捕捉的效果更好,VOF方法只能模拟出气泡变形的大致趋势,各个时刻气泡的高度都要比实验值低。

此外,还可以根据气泡的重心位置和纵横比来比较LS方法和VOF方法的模拟效果。

图7 气泡重心位置随时间的变化

图8 气泡纵横比随时间的变化

由图7可以看出,在气泡产生到发生一段变形时,LS方法和VOF方法都可以很好地预测气泡的重心位置,但随着气泡的继续演化,VOF模拟得到的气泡重心会不断的偏离实验值,而LS方法模拟得到的重心轨迹与实验吻合较好。

气泡的纵横比是指气泡的最大高度与最大宽度的比值,由图7可以看出,在VOF方法中模拟得到的纵横比会沿着实验值震荡,表现为模拟得到的气泡会在轴向上产生周期性的膨胀和收缩,而这一现象是实验观察中不存在的,而LS方法可以很有效地捕捉气泡演化时的形状以及气泡分离的时间。


更多VirtualFlow仿真案例

以下为更多VirtualFlow软件Level Set 界面追踪方法模拟案例:

图9 液体碰撞壁面反弹过程(LS方法)

图10 壁面膜态沸腾LS数值模拟

图11 膜态沸腾(3D)  

   

来源:多相流在线
OpenFOAM碰撞多相流通用ANSAUM控制曲面VirtualFlow
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-06-29
最近编辑:4月前
积鼎科技
联系我们13162025768
获赞 108粉丝 110文章 300课程 0
点赞
收藏
作者推荐

CES:GPU加速的复杂几何气固流动CFD-DEM建模及其应用

计算流体力学和离散元模型的耦合是模拟稠密颗粒系统的强大工具,然而传统的CFD-DEM在处理具有大量颗粒和复杂几何形状的系统时存在限制。本文开发了一种新型的基于GPU的CFD-DEM模型,用于模拟具有大量颗粒和复杂几何形状的气固流动。提出了一种CFD求解器和DEM求解器之间的新型耦合策略,该策略具有高效性和稳定性。通过与实验测量结果的对比,验证了模型的准确性,并将其效率与以往的CFD-DEM模拟进行了比较。然后,作为示例,该模型被用于模拟考虑复杂风口形状和具有大量颗粒的炼铁高炉回旋区中气固流动的动态行为。SIMPOP一、引言由固体颗粒与其周围流体之间强相互作用所构成的稠密颗粒流体系统,在化学工程、矿物加工、药品生产和能源系统等各种工业过程中发挥着至关重要的作用。实际上,这些系统往往具有大量颗粒和复杂几何形状,如循环流化床、旋风分离器和高炉。在工业规模上理解气固系统流动机制对于过程优化至关重要。因此,研究具有大规模尺寸和复杂边界的气固系统具有重要意义。炼铁高炉中的回旋区是稠密气固流动的经典案例。如图1所示,气体通过风口引入,与焦炭颗粒相互作用,推动颗粒向前和向上运动,形成回旋区。因此,回旋区的特点是与风口相连的空隙边界,其中具有大量颗粒。在高炉实际操作中,这一区域对于气体流动和堆积的焦炭颗粒之间的动量、能量和物质交换具有重要意义,并在确定高炉稳定性方面起着关键作用。因此,全面了解高炉内的回旋区行为至关重要。然而,由于回旋区的非透明性、恶劣条件和高昂成本,要深入了解回旋区已被证明是具有挑战性的。图1.炼铁高炉示意图(左),回旋区(右上)和风口细节(右下)以往的研究表明,可以使用CFD-DEM来模拟气固系统的动态行为,但这些研究局限于小尺度系统,颗粒数量有限且几何形状简化。因此,要准确模拟具有大尺寸和复杂边界的气固系统动力学,需要解决两个问题:大量颗粒的处理和具有复杂几何形状系统的建模,特别是对于回旋区系统。在CFD-DEM中,DEM模拟的计算成本很高,特别是在处理大规模场景时。为了解决这个问题,已经应用了几种模型来提高DEM模拟的效率,如粗粒化方法和多相质点网格模型。CGM采用包裹概念来减少颗粒数量,同时通过解析颗粒间碰撞来确保数值准确性。MP-PIC利用包裹来表示具有相同属性的多个真实颗粒的集合,通过引入固相应力模型来简化碰撞。然而,CGM会导致在颗粒尺度上详细的局部信息丢失,以及对于复杂几何形状的流体上局部结构表示的缺失。同时,MP-PIC由于忽略了颗粒间碰撞而牺牲了数值准确性。当这些方法应用于具有复杂系统的气固流动时,往往会遇到挑战。在模型改进方面,一些基于CPU的并行计算技术被用于加速DEM模拟,如MPI(消息传递接口)和OpenMP,以及MPI和OpenMP的耦合。由于进程/核心之间的通信和同步限制,使用OpenMP和MPI的DEM加速比小于线程/核心数,限制了其在大规模DEM模拟中的应用。然而,通过图形处理单元可以显著提高DEM模拟的计算性能。Radeke等首次在GPU上进行了数百万个颗粒的DEM模拟。随后,基于GPU的DEM模拟已应用于许多工业过程,如输送、转鼓、喷动床和流化床。已经证明,基于GPU的DEM在处理大规模计算时表现出很高的效率。第二个问题为模拟具有复杂几何形状的气固流动。目前,大多数CFD-DEM应用都局限于相对简单的系统,因为可以使用CFD中的结构化网格,并在颗粒-壁面碰撞检测中获得更高的计算效率。然而,实际工业中的气固系统通常具有复杂边界。因此,已经开发了一些模型来检测颗粒与复杂几何形状之间的相互作用。例如,Hu等开发了一种颗粒与三角网格边界之间的接触检测模型。Su等提出了一种接触检测算法,可以检测具有多边形网格边界的复杂几何表面上的颗粒接触状态,并进一步提高了计算效率。Chen等也进行了类似的工作。而这些算法都是基于CPU计算实现的,对于大规模计算还需要进一步改进。这在模拟高炉回旋区时尤为重要,因为回旋区具有包括风口在内的复杂边界(如图1c所示)。然而,以往对于高炉回旋区的CFD-DEM模拟并没有考虑到风口的详细特征,通常采用简化的表示方式,如圆面入口、圆柱入口或矩形入口。这些简化的风口模型偏离了工业中实际使用的风口配置,而忽略风口的特定细节可能会导致误导性的模拟结果。因此,在未来的模拟中考虑真实的风口几何形状至关重要,从而实现准确捕捉高炉中气固流动的行为。在这项工作中,开发了一种高效的CFD-DEM方法,通过结合任意边界壁面和大颗粒数来模拟具有大尺度尺寸和复杂边界的气固系统。在该模型中,颗粒模拟通过独立的内部DEM求解器在GPU上执行,并通过引入一种稳定的耦合方法耦合到商业CFD软件ANSYSFluent。随后,对开发的模型进行了验证,并与先前的实验和模拟进行了比较,然后应用于研究考虑大量颗粒和详细风口几何形状的高炉回旋区动力学。二、数值方法2.1.模型描述基于完善的牛顿运动定律,使用离散元法模拟颗粒运动。颗粒i的平移和旋转运动的控制方程给出:其中,和分别是颗粒i的平移速度、角速度、惯性矩和质量。是从颗粒i的中心指向接触点的半径向量。是滚动摩擦系数。是颗粒i的单位角速度。和分别是颗粒i与其邻近颗粒j之间的法向接触力和切向接触力。这些方程的详细信息见表1(Zhu等,2022)。其中:和分别是颗粒体积和直径,是颗粒滑动摩擦系数,和分别是法向和切向重叠量,是颗粒的杨氏模量,是泊松比。和分别是流体粘度和密度,是流体相的体积分数,是流体速度。表1.模拟中使用的力模型基于欧拉模型对流体相进行模拟。其质量和动量守恒的控制方程为:其中,和分别是流体密度、平均速度和脉动速度、压力、流体剪切应力张量和重力加速度。是流体相的体积分数。在模拟中,雷诺应力张量项通过k-ε模型求解。是在给定位置和时间上,所有颗粒对流体施加的每单位体积的合力,由下式给出:其中,是流体单元的体积,N是流体单元中的总颗粒数。2.2.模型实现在这项工作中,模拟是基于内部DEM代码和ANSYS-Fluent的框架进行的,其中颗粒运动通过内部DEM代码在GPU上计算,流体流动则由CPU上的ANSYSFluent解决。CFD与DEM之间的耦合计算(如确定颗粒的CFD单元位置、考虑其体积分数以及将相应的源项合并到动量方程)基于DDPM方案中的函数。DDPM通过将离散相的体积分数引入质量和动量守恒方程中,解决了标准DPM的局限性,允许在稠密颗粒系统中考虑颗粒与连续相之间的相互作用。两个求解器之间的通信通过WindowsWinsock编程实现。1)CFD在CPU上的实现CFD求解器利用MPI在CPU上进行并行化,将计算域划分为多个子域,并将每个子域分配给不同的线程。每个核心都与其他线程同时在其各自的数据集上执行相同的程序。具体来说,在ANSYSFluent中,进程由Cortex、主机进程和一组n个线程的进程组成。主机进程接收来自Cortex的命令,通过process-0向所有进程发送命令,并从process-0收集消息并执行写入文件、数据传输和在所有数据上显示消息等操作。CFD中的计算流程如图2所示。在模拟中,颗粒位置/属性和套接字初始化首先由DEFINE_INIT的UDF(用户定义函数)在主机进程中读取并运行,并且流体流动更新在线程进程中执行。之后,根据从DEM获得的颗粒信息,在主机进程中准备DPM写入并将其分发到所有线程。随后,在线程进程中执行DDPM方案,包括确定颗粒的CFD单元位置,考虑其体积分数,并将相应的源项合并到动量方程。之后,存储在不同线程进程中的颗粒阻力通过node-0传输到主机进程,然后将颗粒信息发送到主机进程中的DEM。最后,在所有进程中清除DPM写入,以避免CFD计算时的中重复计算,为下一次循环迭代做好准备。在这一部分中,使用了UDFDLLEXP的两个函数,“fl_dpm_pre_solve_update”函数用于清除DPM写入,而“fl_dpm_post_solve_update”函数用于创建写入、发送/接收数据和计算阻力。需要注意的是,通过利用ANSYSFluent中内置的DDPM(密集离散相模型)方案,模拟可以有效利用结构化和非结构化网格。与通过UDFDEFINE_SOURCE直接添加动量源项相比,该方法在CFD模拟中提供了更好的收敛性和稳定性,特别是对于涉及复杂几何形状和非结构化网格的情况。主要是由于DDPM方案中使用了基于节点的平滑算法,可以通过平滑加载将颗粒的影响均匀分布到相邻网格节点和相邻单元上,从而减少了网格依赖性和计算不稳定性。图2.CFD-DEM耦合计算流程,方框内的步骤在GPU上执行2)DEM在GPU上的实现图2还展示了DEM在GPU上的计算流程。计算从数据初始化开始,从输入文件读取颗粒属性,然后进行套接字初始化并建立与CFD中套接字的连接。随后,根据输入数据在CPU和GPU中创建阵列的内存,并根据颗粒大小生成颗粒搜索网格。然后,DEM中使用的边界壁面被离散化并分配给颗粒搜索网格。需要注意的是,上述模拟处理在模拟期间仅执行一次。之后,一旦从CFD接收到颗粒属性,就会将其复制到GPU内存,然后在GPU上执行DEM计算,其中包括基于颗粒搜索网格创建颗粒邻居列表,计算颗粒之间以及颗粒与壁面之间的相互作用,以及更新颗粒的位置和速度。DEM计算完成后,将DEM时间与CFD时间进行比较。如果二者同步,颗粒信息从GPU复制到CPU并发送到CFD。否则,将进入新的DEM循环,直到与CFD时间步同步。·将离散化的边界壁面分配给DEM单元在GPUDEM中,一系列三角面元用于表示具有复杂几何的边界壁面以进行高效计算,每个三角面元都被分配一个唯一的编号并与DEM单元相关联。如图3所示,如果三角面元的离散点落在相应的DEM单元内,则将该三角面元添加到该单元的数组中。离散点之间的最小距离设置为DEM单元大小的1/10和最大面长度的最小值,以确保每个三角面元都被正确分配给DEM单元。需要注意的是,该算法仅在主模拟循环之前执行一次,并将排序后的数组存储在GPU内存中以便后续使用。三角面元分配的伪代码可以在补充材料的表S1中看到。图3.DEM单元中的边界壁面和壁面的离散化,用于识别所占用的壁面的DEM单元。·GPU上颗粒与边界壁面的受力计算颗粒与三角面元所组成的边界壁面之间力的计算遵循以下步骤。首先,根据颗粒的位置确定颗粒所在DEM单元的索引。然后,在26个相邻DEM单元以及该单元内包含的所有三角面元上执行循环。在循环过程中,检查颗粒与三角面元之间的距离,并确定颗粒与壁面之间的接触类型。最后,根据接触类型计算颗粒-壁面之间的作用力。关于颗粒与壁面之间接触检测的详细信息可以在Su等的工作中看到。颗粒-壁面作用力计算的伪代码可以在补充材料的表S2中看到。三、模型验证和效率验证3.1.流动动力学验证根据实验室规模二维填充床实验测量中的空腔形态进行验证。填充床的尺寸为长260毫米、高600毫米、宽17.7毫米。用于向填充床供应气体的风口直径和深度分别为10毫米和30毫米。床中的颗粒直径为3-4毫米,密度为1,000kg/m³。填充床的高度为350毫米。实验中注入床中的气体密度为1.23kg/m³,粘度为1.78×10^(-5)kg/(ms),流量为16m³/h。在模拟中,使用相似原理对几何尺寸和颗粒直径进行了缩放。缩放后,几何长、高、宽分别为2,600毫米、3,500毫米和88.6毫米。模拟中注入床中的气体密度为0.716kg/m³,速度为232m/s,颗粒直径为40毫米,填充床中的颗粒总数为11,000。颗粒的杨氏模量为2×10^7Pa,颗粒之间的滚动摩擦系数为0.01。恢复系数设置为0.6,颗粒之间的摩擦系数设置为0.63,而颗粒与壁面之间的摩擦系数设置为0.56。颗粒和气体的其他属性与实验保持一致。DEM网格尺寸设置为颗粒直径的2.5倍,相邻颗粒搜索间隙为颗粒直径的0.8倍,并在后续研究中保持一致。CFD网格尺寸为44.3毫米,CFD生成的网格数为7,640。CFD的时间步长为2×10-4秒,DEM的时间步长为1×10-5秒。使用4个CPU核心进行CFD和1个GPU进行DEM的模拟需要20分钟来计算1秒的物理时间。需要注意的是,在以下工作中,模拟在搭载IntelCorei9-13900KCPU和NVIDIAGeForceRTX3090TiGPU的计算机上执行。所有模拟都在两个求解器中都使用了双精度浮点运算。图4显示了模拟和实验测量之间的空腔形态的对比。结果表明,实验和模拟中的空腔形状具有可比性。具体来说,在实验中,空腔的深度为105.1毫米,高度为120.0毫米,高度与深度的比值为1.14。相比之下,在模拟中,空腔的深度为0.97米,高度为1.196米,高度与深度的比值为1.23。实验与模拟之间的差异约为8.1%,表明模型在空腔形态方面得到了验证。图4.相同条件下模拟与实验之间空腔形态的对比(a)实验(Lu等,2020),(b)模拟和(c)模拟中心面的体积分数3.2.GPUCFD-DEM效率验证为了验证使用基于GPU的DEM-CFD模拟的加速效果,将模拟结果与之前使用CPU进行的模拟进行了比较。该模拟基于实验室规模的鼓泡流化床实验,模拟条件的详细信息如表2和图5所示。模拟和实验都使用了总质量为1.9千克的尼龙珠。表2.性能测试的模拟参数图5.(a)用于性能测试的几何形状和(b)网格图6显示了在气体速度为2.19m/s时颗粒的平均水平/垂直速度。测量区域高度为0.0672m,横向均匀分布5个监测点。每个区域的高度、长度和宽度分别为0.0457m、0.0457m和0.003m。平均速度采用物理时间5s后的平均值,以减少启动效应的影响。图6a-b显示,在相应区域,无论是水平速度还是垂直速度,都与实验数据相当。此外,模拟中床层在0.0413m至0.3461m高度之间的压降为0.70kPa,这也与实验结果的0.69kPa相当。图6.模拟结果与实验数据(Gopalan等,2016)的比较气体速度为2.19m/s时的(a)时间平均颗粒水平速度和(b)时间平均颗粒垂直速度为了验证基于GPU的CFD-DEM模拟的效率,与之前使用基于CPU并行化的粗粒度方法的模拟进行了比较进行了比较。该模拟是在NCI的8核超级计算机Gadi上使用MFIX-DEM软件进行的。图7展示了不同方法在物理时间为10秒时的计算时间。结果表明,基于GPU的模拟在计算时间上大约比基于CPU的模拟低8.5倍。此外,尽管粗粒度方法使用了11,618个颗粒来模拟包含92,948个颗粒的实验室规模的鼓泡流化床,但基于GPU的模拟在计算时间上甚至超过了粗粒度比为2的粗粒度方法。需要注意的是,目前的模拟使用的是基于游戏的GPU,该GPU并未针对双精度计算进行优化。如果使用专业显卡,计算效率将大大提高。图7.不同方法的物理时间为10秒时的计算时间基于CPU的CFD-DEM和粗粒度方法(CGM)的计算时间,数据来自参考文献Wang和Shen(2022)。四、模型应用:高炉回旋区中的气固流动本节通过模拟高炉回旋区来展示模型在处理具有大量颗粒和复杂几何形状的风口的有效性和能力。分析了风口角度对回旋区动力学的影响,并考虑了风口的详细情况,预测了基于颗粒碰撞的风口磨损。4.1.模拟条件图8显示了模拟中使用的几何形状和网格。整个高炉(下部)由40个风口组成。为了降低计算成本,仅模拟了高炉的一部分,代表整体的1/40,并且包含单个风口。此外,模拟排除了受气流影响较小的死区,只关注活性焦炭区域。考虑了三个具有不同倾斜角度的风口。向下倾斜5°(即-5°风口角度)的风口几何细节与工业中使用的相似。水平风口和10°风口的长度与-5°风口相同,而其他尺寸和几何特征则根据-5°风口进行调整。模拟中使用的尺寸和几何细节的详细信息如图8b所示。图8.(a)高炉整体部分和不同角度的风口(b)几何细节(槽型高炉),(c)CFD网格和(d)模拟中使用的初始颗粒图8c展示了模拟中使用的网格。总共使用了9,100个多面体网格,其中风口附近网格尺寸为0.05米,其他区域为0.08米。模拟中共使用了550,000个直径为20毫米、密度为850kg/m^3的焦炭颗粒。在计算之前,所有颗粒都在高炉中随机生成,然后在重力的影响下降落,形成了高度为1,900毫米的填充床,直到系统的动能耗散至接近零,如图8d所示。最后,气体通过风口入口引入高炉,形成回旋区。高炉模拟的参数细节如表3所示。表3.高炉模拟参数4.2.模型的有效性图9展示了总物理时间为5秒的高炉回旋区CFD-DEM模拟中计算时间的百分比。模拟在IntelCorei9-13900K的8核CPU和NVIDIAGeForceRTX3090Ti的GPU上进行。5秒模拟的总计算时间为12.3小时。值得注意的是,CFD求解器占据了大部分计算时间,占总计算时间的84%,而DEM求解器仅占16%,表明CFD求解器是CFD-DEM模拟中的主要瓶颈。此外,在DEM求解器中,创建颗粒邻居列表和计算颗粒-颗粒相互作用占据了大部分计算时间。相比之下,颗粒-壁面相互作用所花费的计算时间不到20%,显示出基于GPU的颗粒-壁面接触检测算法的高效性。图9.不同求解器和关键DEM子程序的计算时间百分比。4.3.回旋区动力学图10显示了不同风口角度的回旋区动力学,其中空隙率等值面为0.55,已证明空隙率为0.55可以准确捕捉回旋区空腔。结果表明,不同风口角度下回旋区的形成是相似的。最初,回旋区沿着垂直于风口的方向扩展,随后逐渐深入高炉内部,导致回旋区高度降低。最终,焦炭颗粒和气体之间实现动态平衡,导致所有风口角度的回旋区在0.5秒后保持稳定。有趣的是,回旋区表现出与风口尖端相对应的更大对齐度,并具有与风口尖端相同的倾斜角度。此外,回旋区高度随着风口角度的增加而减小。图10.不同风口角度(a)水平,(b)-5°和(c)-10°下回旋区的形成。回旋区形状是使用等值面方法生成的,空隙率为0.55。已知回旋区的尺寸对高炉的效率和生产率有显著影响。为了评估不同风口角度下回旋区的尺寸,图11比较了不同风口角度下回旋区尺寸随时间的变化。图11a展示了不同风口角度下回旋区深度随时间的变化,回旋区深度最初迅速增加,然后在0.5秒后趋于稳定,观察到最大回旋区深度出现在-5°风口,但仅略大于-10°的风口得到的结果。图11b-e显示了不同风口角度下高炉内回旋区高度、宽度、表面积和体积的演变,最初都表现出增加,随后随时间减少。此后,略有增加,最终在0.5秒后达到稳定状态。然而,不同风口角度下回旋区尺寸的变化趋势存在显著差异。具体来说,随着风口角度的增加,回旋区高度减小,而回旋区宽度在不同风口角度下趋于相似。此外,在-5°风口角度下观察到最大的表面积,但仅略大于水平风口下所得到的表面积。类似地,回旋区在-5°的风口处达到其最大体积,但回旋区体积明显大于水平风口。此外,-5°风口的回旋区尺寸在0.5秒后的稳定阶段显示出较小的波动,表明与其他两个风口角度相比,回旋区形成更为稳定。上述结果表明,高炉中-5°风口相比其他风口配置具有更好的性能。图11.在空气体积分数为0.55的等值面下,不同风口角度下回旋区特性的演变:(a)深度,(b)高度,(c)宽度,(d)表面积和(e)体积图12显示了在1秒到5秒的回旋区稳定期间的平均尺寸。结果表明,无论是回旋区深度还是宽度,都在-5°风口处达到最大值,并且在-10°风口也观察到相似的值。同时,随着风口角度的增加,回旋区高度逐渐减小。此外,与图11中的结果一致,-5°风口处获得了最大的回旋区体积,而进一步增加风口角度会导致回旋区体积减小。图12.不同风口角度下的平均回旋区尺寸4.4.流场和颗粒动力学图13显示了不同风口角度下的中心截面和水平截面的风速。结果表明,随着风口角度的增加,高风速区域扩大。此外,在水平风口顶部存在一个停滞区域,定义为低速区域(在模拟中小于1m/s)。该停滞区域随着风口角度的增加而减小,并且在风口下方的区域也观察到类似的趋势,这表明较大的风口角度增强了高炉内的流动性。图13.(a)水平风口、(b)-5°风口和(c)-10°风口的中心截面(上)和水平截面(下)的风速分布。图14显示了具有三个颗粒厚度的中心平面上的颗粒速度分布。结果表明,对于不同的风口角度,回旋区中心的颗粒速度较高,而回旋区边界附近的速度相对较低。该现象可能是由于颗粒频繁碰撞而导致的在向前移动时速度下降。然而,随着风口角度的增加,高速颗粒可以达到更低的位置。图14还突出显示了风口附近的颗粒速度矢量,表明回旋区有两个循环——一个在风口的顶部,另一个在风口的底部—导致颗粒不断被送入回旋区。图14.(a)水平风口、(b)-5°风口和(c)-10°风口的中心截面的颗粒速度分布(上)和风口附近的颗粒局部速度矢量(下)。4.5.颗粒活跃区如图14所示,尽管入口速度为173m/s,但大多数颗粒几乎是静止的。为了定义颗粒的活跃区,基于颗粒参与回旋区所需的最小速度,设定了一个0.05m/s的颗粒速度阈值(图14)。图15a展示了在物理时间从1秒到5秒范围内,不同风口角度下的颗粒活跃区。结果表明,颗粒活跃区是一个近球形的形状,其高度比回旋区更高,对于不同的风口角度,都近似于一个椭球体。此外,风口尖端附近的活跃区对于不同的风口角度都倾向于与风口尖端垂直,与回旋区相似。该结果表明风口角度与BF内颗粒活跃区的形状/位置之间存在相关性。随着风口角度的增加,颗粒活跃区的高度降低,说明由于风口倾斜,入口空气提供的能量被消耗在风口尖端底部的颗粒碰撞中。图15b展示了颗粒活跃区的体积和表面积。与回旋区类似,颗粒活跃区在-5°风口处具有最大的体积和表面积,而水平风口的颗粒活跃区体积和表面积最小,这表明风口倾斜增强了BF中的颗粒运动。需要注意的是,如图S1和图S2所示的补充材料中,无论颗粒速度阈值如何变化,颗粒活跃区随风口角度变化的趋势都保持一致。尽管由于阈值定义的不同,形状、体积和表面积会有所变化,但观察到的趋势与图15中所示的趋势相似。图15.在5.0秒时,不同风口角度下(a)颗粒速度为0.05m/s的等值面,(b)颗粒活跃区的平均体积和表面积。4.6.磨损预测已知风口的磨损是高炉风口失效的主要原因之一。为了预测风口的磨损状况,这里采用了广泛使用的Archard磨损模型来测量风口的磨粒磨损,其表达式为:其中,dw和A是几何形状的磨损深度和单元面积,Fn是法向接触力,ΔU是接触时颗粒滑动距离,K是磨损系数,Hv是维氏硬度。通常,风口由铜制成,铜的磨损系数通常为6.0×10^-4,维氏硬度为3.5×10^8Pa。图16a显示了风口上的局部磨损深度;可以看出,对于不同的风口角度,磨损发生在风口末端(风口尖端)附近,但具体的磨损模式因风口角度而异。在水平风口的情况下,最严重的磨损位于风口尖端的大部分区域,而风口尖端底部的磨损相对较低。对于-5°风口,严重的磨损出现在风口尖端的底部和顶部的小部分区域。模拟结果与行业报告中风口失效模式的描述相类似,该报告指出风口尖端顶部和底部的损坏最为严重,如图16b和图16c所示。类似的趋势也在-10°的风口中发现。结果表明,随着风口角度的增加,磨损减少。水平风口在风口尖端顶部和底部区域的显著磨损归因于这些区域存在循环区和风口尖端与颗粒之间的碰撞角度,如图14所示。风口尖端顶部和底部区域较高的碰撞角度导致这些区域磨损严重。通过采用倾斜风口,可以减少颗粒与风口壁之间的碰撞角度,从而随着风口角度的增加而减少磨损。图16.(a)5秒时不同风口角度下风口的局部磨损深度预测,(b)工业中风口尖端的局部缺陷(Rittenschoberetal.,2015),其中颜色越深表示风口损坏的频率越高,(c)工业中风口尖端不同位置的损坏百分比(Rittenschoberetal.,2015),以及(d)风口末端随时间变化的平均磨损深度,线条是对应数据的线性拟合。平均磨损深度是风口尖端区域的平均值。图16d显示了平均磨损深度的定量分析。它显示了在0.5秒之前磨损深度迅速增加,之后对于不同的风口角度,磨损深度随时间线性增加。水平风口、-5°风口和-10°风口的斜率分别为3.88×10-9、3.77×10-9和3.62×10-9,所有风口角度的截距均为2.4×10-9。随着风口角度的增加,斜率减小,表明在相同入口速度下,随着风口角度的增加,磨损减少。值得注意的是,尽管在5秒时磨损深度看似最小,但考虑到风口的寿命,可能会变得相对显著。通常,风口在高炉中的平均寿命为281天,风口尖端处水平、-5°和-10°风口的平均磨损深度分别为11.4毫米、11.0毫米和10.6毫米。此外,通过计算在5秒内的最大磨损并假设这种最大磨损随着时间线性增加,水平、-5°和-10°风口在281天后的最大磨损深度分别高于68.4毫米、66.0毫米和63.5毫米,表明风口损坏主要是由颗粒对风口的磨损造成的。五、结论在本研究中,开发了一种基于CPU-GPU跨平台的CFD-DEM耦合方法,用于模拟具有复杂几何形状的大规模气固流动。颗粒运动由基于GPU的DEM求解器模拟,而流场则在CPU上使用ANSYSFluent进行计算;CFD和DEM之间的耦合通过Winsock数据传输实现。然后应用该模型在颗粒尺度上模拟高炉回旋区动力学。可以得出以下结论:(1)该模型在空腔形态方面得到了验证。该模型的有效性在小型流化床中得到了验证,并且计算时间与CGM方法使用的先前模拟进行了比较,显示出与实验结果良好的一致性以及与CGM方法相比优越的计算速度。(2)CFD计算时间占CFD-DEM模拟的84%,而DEM计算仅占16%。在GPU上的DEM模拟中,大部分时间用于计算颗粒-颗粒力和创建邻居列表。颗粒-壁面力的计算仅占DEM计算时间的16%。(3)在回旋区形成过程中,回旋区的形状与风口尖端对齐,并且具有与风口尖端相同的倾斜角度。此外,当风口向下倾斜5°时,回旋区显示出其最大体积和颗粒活跃区。此外,使用倾斜风口的高炉中,气体和颗粒的流动性都得到了增强。(4)对于水平风口,磨损累积在顶部风口尖端;当风口倾斜时,磨损移动到底部。值得注意的是,最严重的磨损发生在倾斜风口的风口尖端底部,与工业风口中观察到的故障位置一致。随着风口角度的增加,磨损减少,表明较大的风口角度可以更好地防止颗粒碰撞。总体而言,开发的CFD-DEM模型证明了其在模拟复杂几何形状的稠密气固流动的有效性。该模型在研究高炉回旋区动力学方面的成功应用凸显了其在以大量颗粒和复杂几何形状为特征的其他工业过程中的应用潜力。翻译转载自:《ChemicalEngineeringScience》“GPU-acceleratedCFD-DEMmodelingofgas-solidflowwithcomplexgeometryandanapplicationtoracewaydynamicsinindustry-scaleblastfurnaces”来源:多相流在线

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈