CN109522589B - 用于模拟管道颗粒两相流动的非解析方法及电子设备 - Google Patents

Abstract

本公开提供了一种用于模拟管道颗粒两相流动的非解析方法及电子设备。其中所述方法包括：遍历颗粒周边的网格；计算颗粒附近流场的平均孔隙率与速度；计算颗粒所受的拖曳力大小；进行颗粒与流场的耦合，得到颗粒的运动情况和流场的流动情况。本公开对一种用于模拟管道颗粒两相流动的高精度非解析算法能在颗粒大小与网格尺寸相近的情况下很好地描述流场中的颗粒运动，降低了CFD‑DEM耦合算法对网格尺寸的要求，可用于处理液固、气固流化床中的两相流动问题。

Description

用于模拟管道颗粒两相流动的非解析方法及电子设备

技术领域

本公开属于计算力学技术领域，更具体的涉及一种用于模拟管道颗粒两相流动的非解析方法及电子设备。

背景技术

流固两相流动广泛存在与自然环境与工业生产中，例如河流中的沙石沉降、化工领域的喷动床等。采用计算力学的手段获得流场和颗粒的运动情况具有很高的实际应用价值，计算结果可用于预测自然形貌的变化、指导工业设备的设计等。计算流体动力学-离散元(CFD-DEM)耦合算法是一类用于计算流固两相流动的常用方法，采用纳维尔-斯托克斯(Navier-Stokes)方程计算流场的流动，采用离散元法计算颗粒的运动，并通过耦合算法得到流体与固体之间的相互作用。CFD-DEM耦合算法又可细分为两类：解析的耦合算法与非解析的耦合算法。在解析的耦合算法下，流体网格的尺寸需要远小于颗粒的大小，并且在颗粒表面附近需要进行动态加密，计算量较大，无法用于计算成千上万个颗粒的运动情况。而在非解析的耦合算法下，网格尺寸一般是颗粒直径的三倍以上，用拖曳力模型表征颗粒与流场之间的相互作用，计算量远小于解析的耦合算法，所以常用非解析的CFD-DEM耦合算法计算多个颗粒在流场中的运动。

然而，传统的CFD-DEM非解析耦合算法要求网格尺寸在颗粒直径的三倍以上，但是大尺寸的网格很难描述流场中的流动细节，而这些细节信息的缺失不仅会影响我们对流动情况的认知，而且可能使颗粒运动的计算结果与实际情况有所偏差。另外，在计算喷动床、流化床等管道内的流固两相流动时，如果管道内径与颗粒尺寸在同一数量级，计算流体的流动也需要一定数量的网格，那么划分的网格尺寸很难达到颗粒直径的三倍以上。在网格尺寸与颗粒直径相近的条件下，传统的CFD-DEM非解析耦合算法很难准确预测计算区域内颗粒的运动情况。

发明内容

(一)要解决的技术问题

基于以上问题，本公开的目的在于提出一种用于模拟管道颗粒两相流动的非解析方法和电子设备，用于将CFD-DEM非解析耦合算法的适用范围拓展至网格尺寸与颗粒直径相近的情况。

(二)技术方案

为了达到上述目的，作为本公开的一个方面，提供了一种用于模拟管道颗粒两相流动的非解析方法，包括：

设置流场和流场内颗粒的初始条件和边界条件，对流场区域划分网格；

采用希尔伯特曲线遍历颗粒周边的网格；

计算颗粒附近流场的平均孔隙率与速度；

计算颗粒所受的拖曳力大小；

根据所述拖曳力，进行颗粒与流场的耦合，得到颗粒的运动情况和流场的流动情况。

在一些实施例中，用希尔伯特曲线对流场中的网格进行重新编号，并沿该曲线对颗粒周边的网格进行遍历。

在一些实施例中，其中，计算平均孔隙率和速度都采用加权平均的形式，其表达式如下：

式中，

是平均孔隙率，i是网格编号，α_i是第i个网格的孔隙率，W_α是计算平均孔隙率时所用的权函数，r_i是颗粒与第i个网格的中心之间的距离，ΔV_i是第i个网格的体积，

是颗粒附近流场的平均速度，U_f，i是第i 个网格内的流场速度，W_u是计算流场平均速度时所用的权函数。

在一些实施例中，计算平均孔隙率所用的权函数W_α为阶跃函数，其表达式如下：

式中，r是网格中心与颗粒中心之间的距离，L是流场附近球形区域的半径。

在一些实施例中，计算平均速度所用的权函数W_u为高斯函数，其表达式如下：

式中，r是网格中心与颗粒中心之间的距离，σ是用于确定各网格权重的量，与颗粒半径相关。

在一些实施例中，拖曳力采用Abraham公式进行计算，其表达式如下：

式中，C_d是拖曳力系数，Re是雷诺数，F_drag是拖曳力，ρ_f是流体密度，A是颗粒的投影面积，U_r是流场关于颗粒的相对速度大小，d是颗粒的直径。

在一些实施例中，流场采用Navier-Stokes方程进行计算，颗粒运动采用离散元法进行计算，每计算一定时间步长之后，进行两相耦合。

根据本公开的另一方面，提供一种电子设备，包括：

存储器，用于存储可执行指令；以及

处理器，用于执行存储器中存储的可执行指令，以执行如下操作：

设置流场和流场内颗粒的初始条件和边界条件，对所述流场区域划分网格；

采用希尔伯特曲线遍历颗粒周边的网格；

确定颗粒附近流场的平均孔隙率与速度；

确定颗粒所受的拖曳力；

根据上述拖曳力，进行颗粒与流场的耦合，确定颗粒的运动情况和流场的流动情况。

(三)有益效果

(1)本公开用于模拟管道颗粒两相流动的高精度非解析算法可用于计算网格尺寸与颗粒大小相近的情形，所得结果比传统CFD-DEM耦合算法更接近实验结果。

(2)本公开用于模拟管道颗粒两相流动的高精度非解析算法既可以准确地预测单个颗粒在流场中的运动，又可以很好地描述大量颗粒表现出的整体行为。

(3)本公开用于模拟管道颗粒两相流动的高精度非解析算法采用希尔伯特曲线对颗粒周围的网格进行遍历，时间复杂度低，效率高于传统搜索算法。

附图说明

图1为本公开实施例用于模拟管道颗粒两相流动的高精度非解析算法流程图；

图2为本公开实施例1划分的网格以及网格与颗粒相对大小情况；

图3为本公开实施例1所用遍历全场的希尔伯特曲线；

图4为本公开实施例1传统算法与本算法所得颗粒速度与实验结果的比较；

图5为本公开实施例2用本算法得到的两种颗粒在流化床中分层的现象；

图6为本公开实施例2传统算法与本算法所得孔隙率随高度的分布情况与实验结果的比较。

图7为本公开实施例的硬件布置框图。

具体实施方式

为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开作进一步的详细说明。

本公开实施例提供了一种用于模拟管道颗粒两相流动的非解析方法，如图1所示，用于模拟管道颗粒两相流动的非解析方法包括以下步骤：

步骤1、设置流场和流场内颗粒的初始条件和边界条件，对所述流场区域划分网格；具体可以包括：在计算区域内划分用于计算流场的网格，网格尺寸与颗粒大小相近。设置流场与颗粒的物理信息与初始信息，物理信息包括颗粒密度ρ_p、颗粒直径d、流体密度ρ_f、流体运动黏度ν等，初始信息包括颗粒的初始位置(x₀，y₀，z₀)、初始速度U₀等。

步骤2、采用希尔伯特曲线遍历颗粒周边的网格；具体可以包括：生成能够遍历所有网格的希尔伯特曲线，沿希尔伯特曲线对所有网格进行编号。在搜索颗粒周围网格时，沿希尔伯特曲线进行搜索。由于希尔伯特曲线具有分形特性，空间上相邻的网格在希尔伯特曲线上的位置也很接近，沿希尔伯特曲线搜索较少的次数即可遍历颗粒周围的大部分网格。

步骤3、确定颗粒所受的拖曳力大小：可以搜索颗粒周边的网格，计算颗粒所在局部流场内的平均速度和平均孔隙率。利用Abraham公式得到拖曳力F_drag。

然后进行颗粒与流场的耦合，确定颗粒的运动情况和流场的流动情况。具体可以包括以下步骤：

步骤4、用Navier-Stokes方程计算流场，用离散元方法计算颗粒运动。

流体的控制方程为：

式中，α是孔隙率，ρ_f是流体密度，U_f是流体速度，t是时间，p是压强，τ是应力张量，g是重力加速度，∑F_drag是颗粒所受的拖曳力。

颗粒的控制方程为：

式中，m_p是颗粒的质量，C_A是附加质量系数，m_f是与颗粒等体积流体的质量，U_p是颗粒的速度，F_Archimedes是浮力，F_drag是拖曳力，F_p-p是颗粒之间的作用力，F_p-w是颗粒与墙之间的作用力。

对于球形颗粒，上式可写为

式中，d是颗粒的直径，ρ_p是颗粒的密度。

步骤5、每隔一定时间步长，将流场与颗粒耦合，更新颗粒附近流场的平均孔隙率、平均速度，更新颗粒所受的拖曳力，并输出当前时刻的流场信息和颗粒信息。流场信息包括孔隙率α、流体速度U_f等物理量的分布，颗粒信息包括颗粒位置(x，y，z)、颗粒速度U_p等。

下面结合实施例对本公开实施例的方法作详细说明。

实施例1

本实施例以窄管道内的单颗粒沉降来说明本公开用于模拟管道颗粒两相流动的高精度非解析算法。

步骤1、以窄管道内的单颗粒沉降为例。容器内部的长、宽、高为 100mm×100mm×160mm，颗粒直径为15mm，颗粒密度为1120kg/m³，流体为硅油，密度为960kg/m³，动力黏度为0.058N·s/m²，颗粒从容器中轴线上高度为120mm处静止释放。如图2所示，将流场区域划分为 7×7×16的规则六面体网格，每个网格的尺寸与颗粒大小相近。

步骤2、如图3所示，生成遍历所有网格的希尔伯特曲线。

步骤3、确定当前颗粒所在的网格。沿希尔伯特曲线搜索颗粒周边的网格，获得颗粒附近局部流场的信息，计算局部流场的平均孔隙率和平均速度。代入Abraham公式，得到颗粒所受拖曳力。

步骤4、用Navier-Stokes方程计算流场，用离散元方法计算颗粒运动。

步骤5、每隔一定时间步长，将流场与颗粒耦合，更新颗粒附近流场的平均孔隙率、平均速度，更新颗粒所受的拖曳力，并输出当前时刻的颗粒速度U_p。如图4所示，相较于传统算法，本算法所得结果更平稳，与实验结果符合得更好。

实施例2

本实施例以多颗粒流化床来说明本公开用于模拟管道颗粒两相流动的高精度非解析算法。

步骤1、以多颗粒流化床为例。流场区域是直径为104mm、高为200mm 的圆柱体。流体为水，密度为1000kg/m³，运动黏度为1.0×10^-6m²/s。颗粒密度为2960kg/m³，有两种直径：2mm和4mm，两种颗粒的质量比为1∶1。将流场区域划分为21×21×40的规则六边形网格，网格尺寸与大颗粒的大小相当，为达到小颗粒直径的三倍。

步骤2、生成遍历所有网格的希尔伯特曲线。

步骤3、对每个颗粒，确定其当前所在的网格。沿希尔伯特曲线搜索颗粒周边的网格，获得颗粒附近局部流场的信息，计算局部流场的平均孔隙率和平均速度。代入Abraham公式，得到颗粒所受拖曳力。

步骤4、用Navier-Stokes方程计算流场，用离散元方法计算颗粒运动。

步骤5、每隔一定时间步长，将流场与颗粒耦合，更新颗粒附近流场的平均孔隙率、平均速度，更新颗粒所受的拖曳力，并输出当前时刻每个颗粒的位置(x，y，z)。

如图5所示，本算法可以算出两种颗粒明显的分层现象，大颗粒501 在下，小颗粒502在上。根据每个颗粒的位置和直径，计算孔隙率随高度的分布情况。如图6所示，相较于传统算法，本算法所得结果与实验结果更吻合。

图7是示出了根据本公开实施例的示例硬件布置700的框图。电子设备700包括处理器706(例如，微处理器(μP)、数字信号处理器(DSP) 等)。处理器706可以是用于执行本文描述的流程的不同动作的单一处理单元或者是多个处理单元。电子设备700还可以包括用于从其他实体接收信号的输入单元702、以及用于向其他实体提供信号的输出单元704。输入单元702和输出单元704可以被布置为单一实体或者是分离的实体。

此外，电子设备700可以包括具有非易失性或易失性存储器形式的至少一个存储器708，例如是电可擦除可编程只读存储器(EEPROM)、闪存、和/或硬盘驱动器。存储器708包括计算机程序710，该计算机程序710包括代码/计算机可读指令，其在由电子设备700中的处理器706执行时使得电子设备700和/或包括电子设备700在内的系统可以执行例如上面结合图 1所描述的流程及其任何变形。

计算机程序710可被配置为具有例如计算机程序模块710A～710E架构的计算机程序代码。因此，在例如电子设备700的计算机程序中的代码包括：模块710A，用于设置流场和流场内颗粒的初始条件和边界条件，对所述流场区域划分网格。计算机程序中的代码还包括：模块710B，用于采用希尔伯特曲线遍历颗粒周边的网格。计算机程序中的代码还包括：模块710C，用于确定颗粒附近流场的平均孔隙率与速度。计算机程序中的代码还包括：模块710D，用于确定颗粒所受的拖曳力大小。计算机程序中的代码还包括：模块710E，用于进行颗粒与流场的耦合，确定颗粒的运动情况和流场的流动情况。

计算机程序模块实质上可以执行图1中所示出的流程中的各个动作。

尽管上面结合图7所公开的实施例中的代码手段被实现为计算机程序模块，其在处理器706中执行时使得电子设备700执行上面结合图1所描述的动作，然而在备选实施例中，该代码手段中的至少一项可以至少被部分地实现为硬件电路。处理器可以是单个CPU(中央处理单元)，但也可以包括两个或更多个处理单元。例如，处理器可以包括通用微处理器、指令集处理器和/或相关芯片组和/或专用微处理器(例如，专用集成电路 (ASIC))。处理器还可以包括用于缓存用途的板载存储器。计算机程序可以由连接到处理器的计算机程序产品来承载。计算机程序产品可以包括其上存储有计算机程序的计算机可读介质。例如，计算机程序产品可以是闪存、随机存取存储器(RAM)、只读存储器(ROM)、EEPROM，且上述计算机程序模块在备选实施例中可以用UE内的存储器的形式被分布到不同计算机程序产品中。

在以上所提供的实施例中，说明了大量具体细节。然而，能够理解，本公开的实施例可以在没有这些具体细节的情况下实践。在一些实例中，并未详细示出公知的方法、结构和技术，以便不模糊对本说明书的理解。

类似地，应当理解，为了精简本公开并帮助理解各个公开方面中的一个或多个，在上面对本公开的示例性实施例的描述中，本公开的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开解释成反映如下意图：即所要求保护的本公开要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，公开方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本公开的单独实施例。

还需要说明的是，本文可提供包含特定值的参数的示范，但这些参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应值。实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向，并非用来限制本发明的保护范围。此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。并且上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。

以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。

Claims (6)

1.一种用于模拟管道颗粒两相流动的非解析方法，包括：

设置流场和流场内颗粒的初始条件和边界条件，对流场区域划分网格；其中，在计算区域内划分用于计算流场的网格，网格尺寸与颗粒大小相近；

采用希尔伯特曲线遍历颗粒周边的网格；

确定颗粒附近流场的平均孔隙率与速度；

确定颗粒所受的拖曳力；

根据之前所述的拖曳力，进行颗粒与流场的耦合，将拖曳力施加至颗粒及流体上，确定颗粒的运动情况和流场的流动情况；

其中，所述确定颗粒附近流场的平均孔隙率与速度，包括：

采用加权平均的形式计算平均孔隙率和速度，其表达式如下：

式中，

是平均孔隙率，i是网格编号，α_i是第i个网格的孔隙率，W_α是计算平均孔隙率时所用的权函数，r_i是颗粒与第i个网格的中心之间的距离，ΔV_i是第i个网格的体积，

是颗粒附近流场的平均速度，U_f，i是第i个网格内的流场速度，W_u是计算流场平均速度时所用的权函数；

所述拖曳力采用Abraham公式进行计算，其表达式如下：

式中，C_d是拖曳力系数，Re是雷诺数，F_drag是拖曳力，ρ_f是流体密度，A是颗粒的投影面积，U_r是流场关于颗粒的相对速度大小，d是颗粒的直径。

2.根据权利要求1所述的方法，其中，用希尔伯特曲线对流场中的网格进行重新编号，并沿该曲线对颗粒周边的网格进行遍历。

3.根据权利要求1所述的方法，其中，计算平均孔隙率所用的权函数所述W_α为阶跃函数，其表达式如下：

式中，r是网格中心与颗粒中心之间的距离，L是流场附近球形区域的半径。

4.根据权利要求1所述的方法，其中，计算平均速度所用的权函数所述W_u为高斯函数，其表达式如下：

式中，r是网格中心与颗粒中心之间的距离，σ是用于确定各网格权重的量，与颗粒半径相关。

5.根据权利要求1所述的方法，其中，流场采用Navier-Stokes方程进行计算，颗粒运动采用离散元法进行计算，每计算设定时间步长之后，进行两相耦合。

6.一种电子设备，包括：

存储器，用于存储可执行指令；以及

处理器，用于执行存储器中存储的可执行指令，以执行如下操作：

设置流场和流场内颗粒的初始条件和边界条件，对流场区域划分网格；其中，在计算区域内划分用于计算流场的网格，网格尺寸与颗粒大小相近；

采用希尔伯特曲线遍历颗粒周边的网格；

确定颗粒附近流场的平均孔隙率与速度；

确定颗粒所受的拖曳力大小；

根据所述拖曳力，进行颗粒与流场的耦合，确定颗粒的运动情况和流场的流动情况；

其中，所述确定颗粒附近流场的平均孔隙率与速度，包括：

采用加权平均的形式计算平均孔隙率和速度，其表达式如下：

式中，

是平均孔隙率，i是网格编号，α_i是第i个网格的孔隙率，W_α是计算平均孔隙率时所用的权函数，r_i是颗粒与第i个网格的中心之间的距离，ΔV_i是第i个网格的体积，

是颗粒附近流场的平均速度，U_f，i是第i个网格内的流场速度，W_u是计算流场平均速度时所用的权函数；

所述拖曳力采用Abraham公式进行计算，其表达式如下：

式中，C_d是拖曳力系数，Re是雷诺数，F_drag是拖曳力，ρ_f是流体密度，A是颗粒的投影面积，U_r是流场关于颗粒的相对速度大小，d是颗粒的直径。

Images