CN113378445B - 一种基于离散模拟的气液多相系统计算方法及系统 - Google Patents

一种基于离散模拟的气液多相系统计算方法及系统 Download PDF

Info

Publication number
CN113378445B
CN113378445B CN202110504761.0A CN202110504761A CN113378445B CN 113378445 B CN113378445 B CN 113378445B CN 202110504761 A CN202110504761 A CN 202110504761A CN 113378445 B CN113378445 B CN 113378445B
Authority
CN
China
Prior art keywords
particle
particles
gas
speed
liquid
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110504761.0A
Other languages
English (en)
Other versions
CN113378445A (zh
Inventor
陈飞国
葛蔚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Process Engineering of CAS
Original Assignee
Institute of Process Engineering of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Process Engineering of CAS filed Critical Institute of Process Engineering of CAS
Priority to CN202110504761.0A priority Critical patent/CN113378445B/zh
Publication of CN113378445A publication Critical patent/CN113378445A/zh
Application granted granted Critical
Publication of CN113378445B publication Critical patent/CN113378445B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Computing Systems (AREA)
  • Artificial Intelligence (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Fluid Mechanics (AREA)
  • Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Data Mining & Analysis (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于离散模拟的气液多相系统计算方法及系统，所述方法包括：根据模拟设置信息设定待模拟的气液多相系统初始状态；对待模拟的气液多相系统进行空间分解得到多个子区域，对每个子区域按照作用截断距离信息进行划分，得到若干个网格；根据初始状态，将粒子归属到不同的网格中，建立每个粒子与网格的映射关系；采用蛙跳算法同步迭代更新每个粒子的速度与位置，进而更新每个粒子与网格的映射关系；再根据更新后的速度与位置重新计算每个子区域不同粒子间的作用，进而由蛙跳算法继续更新每个粒子的速度；反复迭代更新每个粒子的速度与位置，并输出模拟结果。

Description

一种基于离散模拟的气液多相系统计算方法及系统
技术领域
本发明涉及气液多相系统的离散计算技术领域，特别涉及一种基于离散模拟的气液多相系统计算方法及系统，关于气液界面形变和破碎过程高精度模拟计算。
背景技术
气液多相系统指同时具有气体和液体的系统，广泛存在于自然环境和工业生产过程中。气液多相系统的模拟计算重点需要解决两个方面的问题：表达气液两相运动和追踪气液相界面，模拟计算方法主要可以分为基于连续模拟的计算方法和基于离散模拟的计算方法。
基于连续模拟的计算方法将气液两相处理为可互相渗透的两种连续体，根据气液两相在空间中的体积分率分布情况近似计算获得气液相界面，气液相界面的近似计算方法主要有标记网格法、流体体积法、水平集法和波前追踪法等。有关气液多相系统连续模拟计算和气液相界面追踪方法可参见参考文献1：Guan Heng Yeoh and Jiyuan Tu.Chapter 3-Solution Methods for Multi-Phase Flows.In Computational Techniques forMultiphase Flows(eds.Guan Heng Yeoh and Jiyuan Tu),95-242.Oxford:Butterworth-Heinemann,2010.
基于离散模拟的计算方法将气液两相处理为两个分别带有各相性质的粒子群，粒子的运动和粒子之间的相互作用描述了气液多相系统运动，而气液相界面则有两种不同粒子的位置分布自然获得。不同的粒子间相互作用模型构成了不同离散模拟方法。在光滑粒子动力学中，流体离散成许多细小的微元(称作光滑粒子)，微元上的流体力学特性按照周围粒子利用权函数插值得到，邻近微元之间的作用遵循Navier-Stokes方程。光滑粒子间的作用可表达为粒子受力，如
其中ri、rj是i粒子和j粒子的位置，vi、vj是i粒子和j粒子的速度，μi、μj是i粒子和j粒子的黏度；ni、nj是i粒子和j粒子当地的数密度，Pi、Pj是i粒子和j粒子的压力，W(ri-rj,h)是权函数，h是光滑长度。有关光滑粒子动力学计算可参见参考文献2：Guirong Liu andMoubin Liu.Smoothed particle hydrodynamics:a meshfree particle method.Worldscientific,2003.
拟颗粒模型将软球模型的时驱算法结合到硬球模型中。每一个时间步内，所有拟颗粒按照同一时间步长独立运动；在结束时，若与其他拟颗粒重叠即拟颗粒之间的距离|r1-r2|小于它们的半径之和以及r1-r2和v1-v2的点积为负值，则它们将以光滑硬球方式碰撞，
其中v1’、v2’和v1、v2分别是两粒子碰撞后和碰撞前的速度，m1和m2是它们的质量，e是拟颗粒的回复系数，一般为1。下一个时间步，拟颗粒按照新速度移动到新的位置并以此循环。若同一时间步内一个拟颗粒与多个拟颗粒发生碰撞，拟颗粒间的碰撞将按照设定顺序轮换进行，以确保各向同性和空间均匀性。拟颗粒模型能较好表达表现气体运动性质，并且在计算处理上也具备简单高效的优势。有关拟颗粒模型的具体情况可参见参考文献3：Wei Ge and Jinghai Li.Simulation of particle–fluid systems with macro-scalepseudo-particle modeling.Powder Technology,2003,137(1):99-108.
现有技术中存在多种实现气液多相系统模拟的相关算法，但这些算法中大部分以连续流体方式处理气液两相，或者在离散计算方法中仅以不同作用参数区分气液两相。这些算法并不能够准确表达气液运动特征的差异及气液相界面变形和破碎过程。例如：
1、在连续流体方式的计算方法中，气液相界面通过流体体积法或水平集法等界面追踪算法近似获得，该计算过程繁琐、所需计算量较大；在界面存在较大变形甚至破碎过程中，上述界面追踪方式效率低下或者难以完成。
2、现有离散方式计算气液多相系统的方法采用相同的作用方式处理气液两相，而仅以不同的粒子属性及作用参数区分；通常气体行为主要由分子间的排斥作用决定，而液体分子间的吸引作用对保持液体行为也有重要作用，上述方式难以准确表达气液两相运动特征的差异。
发明内容
本发明的目的是克服现有气液多相系统模拟计算中无法准确表达气液行为特征和无法精确高效追踪气液界面的缺陷，提供一种实现气液多相系统模拟及气液界面形变和破碎过程高精度计算的相关方法。
为了实现上述目的，本发明提出了一种基于离散模拟的气液多相系统计算方法，所述方法包括：
步骤1)根据模拟设置信息设定待模拟的气液多相系统初始状态；所述初始状态包括气体粒子数目、液体粒子数目、每个气体粒子的初始位置信息及初始速度，每个液体粒子的初始位置信息、初始速度及初始作用力；
步骤2)对待模拟的气液多相系统进行空间分解得到多个子区域，对每个子区域按照作用截断距离信息进行划分，得到若干个网格；
步骤3)根据初始状态，将粒子归属到不同的网格中，建立每个粒子与网格的映射关系；
步骤4)采用蛙跳算法同步迭代更新每个粒子的速度与位置，进而更新每个粒子与网格的映射关系；再根据更新后的速度与位置重新计算每个子区域不同粒子间的作用，进而由蛙跳算法继续更新每个粒子的速度；
步骤5)判断是否达到模拟结束条件，判断为是，转至步骤6)，否则转至步骤4)；
步骤6)输出模拟结果。
作为上述方法的一种改进，所述模拟设置信息包括模拟区域信息、气体粒子信息、液体粒子信息、气体定向速度、液体注入位置和液体注入速度。
作为上述方法的一种改进，所述采用蛙跳算法同步迭代更新每个粒子的速度与位置，具体包括：
根据粒子在t时刻的速度v(t)、位置r(t)和作用力f(t)，采用蛙跳方式，将速度初步更新至(t+Δt/2)时刻，得到v(t+Δt/2)，从而将位置更新至(t+Δt)时刻得到r(t+Δt)，具体表示如下：
r(t+Δt)＝r(t)+v(t+Δt/2)Δt
其中，m为粒子质量。
作为上述方法的一种改进，所述更新每个粒子与网格的映射关系；具体包括：
先删除原网格内的映射关系，再依据粒子更新后的位置建立粒子到网格的新映射关系。
作为上述方法的一种改进，所述根据更新后的速度与位置重新计算每个子区域不同粒子间的作用，进而由蛙跳算法继续更新每个粒子的速度；具体包括：
针对每个子区域，读取该子区域以及该子区域相邻子区域内的粒子位置信息，采用shift模式，进行不同子区域间的粒子迁移和边界粒子传递；
读取该子区域的一个粒子的位置信息r(t+Δt)，根据更新后的映射关系确定该粒子所在网格，根据所在网格搜寻邻居网格内的其他粒子；
计算该粒子与其他粒子间距离，当距离在截断距离内则构成一个待计算粒子对；
遍历该粒子的每个待计算粒子对，根据粒子对的粒子类型选择相应的粒子作用方式和作用参数，计算对应的作用力或碰撞后的速度更新值；
由该粒子的每个粒子对的作用得到该粒子在(t+Δt)时刻的合力f(t+Δt)；
根据下式得到该粒子在(t+Δt)时刻的速度v(t+Δt)：
作为上述方法的一种改进，所述遍历该粒子的每个待计算粒子对，根据粒子对的粒子类型选择相应的粒子作用方式和作用参数，计算对应的作用力或碰撞后的速度更新值；具体包括：
对于气体粒子间采用拟颗粒作用方式，两个粒子即时更新速度信息；
对于液体粒子间采用光滑粒子动力学作用方式，获得粒子作用力；
对于气液粒子间采用拟颗粒作用方式，即时更新气体拟颗粒的速度信息，并将碰撞过程发生的液体粒子动量变化量除以持续时间转换为作用力；
所述作用参数包括作用截断距离rc、作用特征长度h、黏度系数μ、表面张力系数σ、参考密度ρ0、参考压强p0以及碰撞回复系数e。
一种基于离散模拟的气液多相系统计算系统，所述系统包括：初始设置模块、子区域划分模块、多个并行的处理模块和输出模块；其中，
所述初始设置模块，用于根据模拟设置信息设定待模拟的气液多相系统初始状态；所述初始状态包括气体粒子数目、液体粒子数目、每个气体粒子的初始位置信息及初始速度，每个液体粒子的初始位置信息、初始速度及初始作用力；
所述子区域划分模块，用于对待模拟的气液多相系统进行空间分解得到多个子区域，并分别配置给对应的处理模块；所述处理模块有多个，为并行处理；
所述处理模块，用于对本子区域按照作用截断距离信息进行划分，得到若干个网格；用于根据初始状态，将粒子归属到不同的网格中，建立本子区域每个粒子与网格的映射关系；还用于采用蛙跳算法同步迭代更新每个粒子的速度与位置，进而更新每个粒子与网格的映射关系，再根据更新后的速度与位置重新计算每个子区域不同粒子间的作用，进而由蛙跳算法继续更新每个粒子的速度；并将每次迭代更新的每个粒子的速度与位置发送至输出模块；
所述输出模块，用于输出模拟结果。
与现有技术相比，本发明的优势在于：
1、本发明建立了气液多相系统的离散粒子计算方法，气体粒子和液体粒子的位置差异自然构成了气液相界面，界面确定方式简单直接，且不引入额外的近似处理算法，在界面捕捉精度和效率上均有改善；
2、拟颗粒模型表达的粒子间相互作用为不可叠加的硬球碰撞，需要即时更新碰撞后粒子信息，可以较好表征气相运动特征；而表征液体的光滑粒子动力学为可叠加的软球作用，以可加和的作用力方式施加于液体粒子上，本发明的方法对气液多相系统中涉及的气体粒子和液体粒子采用了完全不同的作用方式，清晰表达和区分气液两相不同的运动特征。
附图说明
图1是在本发明的一个实施例中所提供的并行计算系统的示意图；
图2是本发明的方法的流程图；
图3是本发明的方法中对并行计算区域进行分解的示意图；
图4是本发明在一个实施例中在邻居进程间所采用的shift模式的示意图；其中图4(a)表示shift模式中各进程间不同方向的传输顺序，图4(b)表示shift模式的第一个方向传输，图4(c)表示shift模式的第二个方向传输；
图5是采用本发明的方法对气液多相系统进行模拟的结果示意图；其中图5(a)表示为t在1时的结果，图5(b)表示t在50时的结果，图5(c)表示t在100时的结果，图5(d)表示t在150时的结果，图5(e)表示t在450时的结果。
具体实施方式
下面结合附图和实施例对本发明的技术方案进行详细的说明。
实施例1
本发明的实施例1提出了一种基于离散模拟的气液多相系统计算方法。具体包括：
步骤1)根据模拟设置信息设定待模拟的气液多相系统初始状态；所述初始状态包括气体粒子数目、液体粒子数目、每个气体粒子的初始位置信息及初始速度，每个液体粒子的初始位置信息、初始速度及初始作用力；
步骤2)对待模拟的气液多相系统进行空间分解得到多个子区域，对每个子区域按照作用截断距离信息进行划分，得到若干个网格；
步骤3)根据初始状态，将粒子归属到不同的网格中，建立每个粒子与网格的映射关系；
步骤4)采用蛙跳算法同步迭代更新每个粒子的速度与位置，进而更新每个粒子与网格的映射关系；再根据更新后的速度与位置重新计算每个子区域不同粒子间的作用，进而由蛙跳算法继续更新每个粒子的速度；
步骤5)判断是否达到模拟结束条件，判断为是，转至步骤6)，否则转至步骤4)；
步骤6)输出模拟结果。
本发明在模拟气液多相系统时，虽然现在单个处理器单元的计算性能相当优异，但为了更好地对气液多相系统中大量粒子进行模拟，本发明采用并行计算系统实现对气液多相系统的模拟计算。所述的并行计算系统由多个计算节点组成，在每个计算节点上有若干处理器单元，不同计算节点间通过网络连接。
在本发明的一个实施例中，提供了如图1所示的并行计算系统，该计算系统包括4个计算节点，每个计算节点安装有2块公司生产的Xeon(R)CPU，由Linux系统构建计算环境，采用InfiniBand高速网络连接各计算节点。
在本实施例中，假设用上述的并行计算系统对一种气液多相系统做离散模拟，所要模拟区域的大小为800×240×120(此处的大小为无量纲值，可以根据实际体系需要转化为带单位的真实值，其他物理量同理可转换成真实值)，初始区域内充满气体粒子(总数Ng＝1.332×107)，液体从顶面一圆孔以一定流量和速度注入区域内，在一侧面区域设置气体粒子定向速度为一设定值，从而对液柱进行横向雾化；模拟过程中，随着液体注入及气体运动，模拟区域内的气体粒子数和液体粒子数会发生变化。
下面参考图2，对如何在图1所示的并行计算系统上实现对气液多相系统的模拟进行说明。
将所要模拟的气液多相系统的上述信息读入到前述的并行计算系统中，在所述的计算机系统中实现对模拟区域的空间分解。从前面的说明中已经知道，本发明为了提高效率，将多个具有处理器单元的计算节点通过网络并行连接。为了提高模拟计算的效率，需要将模拟区域中的不同空间分配给不同的处理器单元，这种分配也被称为任务划分。因此，首先需要对模拟区域做空间分解，即将三维的模拟区域做立体分割，得到相应的子区域。在本实施例中，采用二维分割的方式将所述的模拟区域分割成Nx×Ny(本实施例10×8)个子区域，在各个子区域内再根据作用截断距离信息(作用截断距离是指当两个粒子间距离大于该作用截断距离时，则两粒子之间就被视作不发生作用)划分网格。在划分子区域时，所划分子区域的个数通常与并行计算机系统内的处理器单元数相关，一个处理器单元对应一个子区域。本实施例中，由子区域进一步划分所得到的网格的大小为2.55。
在对模拟区域做空间分解后，并行计算机系统根据模拟区域中的粒子所在的位置，在粒子与粒子所在网格间建立对应关系，即建立粒子到网格的映射。通过映射得到的相关信息，可以知道一个粒子在哪个网格中，也可以知道在一个网格中有哪些粒子。
在为处理器单元分配对应网格的粒子后，还需要将作用方式及作用参数、模拟设置等信息保存到内存中。以本实施例所要模拟的气液多相系统为例，该系统具有气、液两相。在气液两相体系中，粒子间两两作用共有3种作用方式及作用参数。所述作用方式及作用参数具体包括：作用方式1，气体-气体间拟颗粒作用，作用截断距离rc＝0.5，碰撞回复系数e＝1.0；作用方式2，液体-液体间光滑粒子动力学作用，作用截断距离rc＝2.55，作用特征长度h＝0.85，黏度系数μ＝1，表面张力系数σ＝2，参考密度ρ0＝2.30，参考压强p0＝295；作用方式3，气体-液体间碰撞作用,作用截断距离rc＝0.675，碰撞回复系数e＝0.9。
要保存到内存中的模拟设置信息则包括：侧面区域设置气体定向速度Ug＝8，液体注入位置为以顶面(40，60)为中心直径20的圆孔，每隔0.1时间注入速度为5厚度0.5的液柱等。
通过以上操作，为各个处理器单元分配任务并赋予相应的参数后，就可以在各个进程中为对应子区域内的粒子执行相应的计算。在计算过程中，由于所要模拟的气液多相系统中的各个粒子是随着时间而不断运动着的，因此需要对粒子的位置和速度信息做更新，了解粒子当前所在的位置和速度的大小。在更新粒子位置和速度时，可采用蛙跳方式，即：
r(t+Δt)＝r(t)+v(t+Δt/2)Δt
具体而言，在更新粒子位置和速度时，首先将粒子在t时的速度v(t)根据当时的受力f(t)更新到v(t+Δt/2)，从而将其位置r(t)更新到r(t+Δt)；然后按t+Δt时的各粒子位置计算受力f(t+Δt)，把速度由v(t+Δt/2)更新到v(t+Δt)，从而完成粒子位置、速度的同步。
在更新粒子的位置后，需要重新实现粒子到网格的映射。在前文中已经提到，粒子的位置在运动过程中会发生变化，可能会从一个网格到了另一个网格，因此，需要对粒子到网格的映射关系加以更新。在更新过程中，需要先删除原网格内的相应映射，再依据粒子新位置建立粒子到网格的新映射关系。
通过上述的操作，完成了粒子到网格映射的更新，粒子到网格映射的更新主要解决了在同一个子区域内的粒子在不同网格间迁移的问题，但在实际应用中，粒子也有可能从一个子区域移动到另一个子区域，即由一个处理器单元转到另一个处理器单元做相关的粒子计算，因此还需要实现在不同处理器单元间粒子的迁移操作。为了方便理解，以图3为例，对进程间粒子迁移的操作进行说明。在图3中，根据前述任务划分的需求，整个模拟区域[GlobalLeftBottom,GlobalRightTop]被分割成Nx×Ny个子区域，一个进程Pi,j(进程编号为i+j×Nx)负责一个子区域内([LocalLeftBottom,LocalRightTop])的计算。在前面已经提到，子区域可以进一步划分为网格，在图4中，用小方块表示网格。子区域按照由外到内的顺序，还可以分为外边界区、内边界区以及子区域内部，其中，一个子区域的外边界区是相邻子区域的内边界区，类似的，一个子区域的内边界区是相邻子区域的外边界区。例如，图3中进程Pi,j所对应子区域左侧外边界区是其左侧相邻的进程Pi-1,j所对应的子区域的右侧内边界区。如果某个粒子从进程Pi,j所对应子区域左侧内边界区运动到了左侧外边界区，实际上也就运动到了进程Pi-1,j所对应的子区域的右侧内边界区，因此对该粒子的处理由进程Pi,j转移到了进程Pi-1,j。
此外，在计算粒子间相互作用力时，由于子区域的边界区域内的粒子需要其他子区域内的粒子信息才能完成完整的受力计算，因此一个进程还需要得到邻居进程的边界粒子信息，这种现象也被称为边界粒子传递。
上述的进程间粒子迁移与边界粒子传递都可以采用相同的方法实现，因此下面做统一的说明：
步骤A、各个线程各自将位于内边界区域和外边界区域中的各个粒子信息保存到内存缓冲区内；
步骤B、各进程将所述发送缓冲区内的信息传输到目标邻居进程的接收缓冲区，将其他邻居进程传入的信息放置于本进程的接收缓冲区内；
步骤C、解析接收到的内存缓冲区内粒子信息，将跨界粒子添加到本进程的粒子信息列表中，或者替换位于外边界网格的粒子信息。
在上述的步骤B中，为了提高通信效率，可以采用现有技术中的shift模式实现相邻进程间的通信。以图3中的空间二维分割为例，shift模式将通信过程分为两个阶段，分别是x方向的传递和y方向的传递，如图4所示。对于空间二维分割而言，shift模式中一个进程只需与相邻四个进程进行通信，与点对点传输模式相比，shift模式通信次数大大减少。
完成关于进程间粒子迁移和边界粒子传递的相关操作后，就可以对粒子间作用进行计算。其实现步骤如下：
步骤1、读入一个粒子信息，根据粒子位置信息确定粒子所在网格，根据对应网格搜寻邻居网格内的其他粒子；
步骤2、计算本粒子与其他粒子间距离，两者距离在截断距离rc内即构成待计算粒子对；
步骤3、根据待计算粒子对的粒子类型选择相应的粒子作用方式和作用参数，计算待计算粒子对间作用：气体粒子间采用拟颗粒作用方式，两个粒子即时更新速度信息；液体粒子间采用光滑粒子动力学作用获得粒子作用力；气液粒子间采用拟颗粒作用方式即时更新气体拟颗粒的速度信息，并将碰撞过程发生的液体粒子动量变化量除以持续时间转换为作用力；
步骤4、在搜寻完一个邻居网格内的所有其他粒子后轮换邻居网格，重复上述的步骤4-4-2)和4-4-3)，直到执行完与周边所有邻居粒子间的作用计算后将包括速度和作用力大小的粒子信息写回。
以上就是粒子间作用计算的基本实现步骤。
在得到粒子所受到的作用力后，可以根据作用力再次更新粒子的速度，对粒子速度更新的实现可参照前面的说明。然后对气液多相系统的模拟是否结束进行判断，若没有结束，则重新执行上述从更新粒子位置起的各个操作步骤，直到模拟结束后将相关信息从并行计算系统输出。以上是对在并行计算系统上实现气液多相系统模拟的基本流程。
采用本发明的方法可以取得良好的模拟效果。对上述实施例中所模拟的气液多相系统的模拟结果在图5中做了展示。图中仅展示了液体粒子的分布情况，填充在空间内的其余部分为气体。在图5(a)中t＝1时，液柱从顶部以圆孔进入区域内；在图5(b)中t＝50时，液柱持续进入区域内，并被横向气流作用偏离初始方向；而在图5(c)中t＝100时，由于气流对液柱的剪切作用，液柱发生破裂产生不同大小液滴；在图5(d)(t＝150)和图5(e)(t＝450)中，液柱继续进入区域内，气流继续作用于液柱使之发生偏向及破裂，形成周期性的雾化运动。此现象与实际气液多相系统的横流雾化典型过程相符，因此采用上述方法所实现的模拟结果较为真实可信。
实施例2
本发明的实施例2提出了一种基于离散模拟的气液多相系统计算系统，所述系统包括：初始设置模块、子区域划分模块、多个并行的处理模块和输出模块；其中，
所述初始设置模块，用于根据模拟设置信息设定待模拟的气液多相系统初始状态；所述初始状态包括气体粒子数目、液体粒子数目、每个气体粒子的初始位置信息及初始速度，每个液体粒子的初始位置信息、初始速度及初始作用力；
所述子区域划分模块，用于对待模拟的气液多相系统进行空间分解得到多个子区域，并分别配置给对应的处理模块；所述处理模块有多个，为并行处理；
所述处理模块，用于对本子区域按照作用截断距离信息进行划分，得到若干个网格；用于根据初始状态，将粒子归属到不同的网格中，建立本子区域每个粒子与网格的映射关系；还用于采用蛙跳算法同步迭代更新每个粒子的速度与位置，进而更新每个粒子与网格的映射关系，再根据更新后的速度与位置重新计算每个子区域不同粒子间的作用，进而由蛙跳算法继续更新每个粒子的速度；并将每次迭代更新的每个粒子的速度与位置发送至输出模块；
所述输出模块，用于输出模拟结果。
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (5)

1.一种基于离散模拟的气液多相系统计算方法，所述方法包括：
步骤1)根据模拟设置信息设定待模拟的气液多相系统初始状态；所述初始状态包括气体粒子数目、液体粒子数目、每个气体粒子的初始位置信息及初始速度，每个液体粒子的初始位置信息、初始速度及初始作用力；
步骤2)对待模拟的气液多相系统进行空间分解得到多个子区域，对每个子区域按照作用截断距离信息进行划分，得到若干个网格；
步骤3)根据初始状态，将粒子归属到不同的网格中，建立每个粒子与网格的映射关系；
步骤4)采用蛙跳算法同步迭代更新每个粒子的速度与位置，进而更新每个粒子与网格的映射关系；再根据更新后的速度与位置重新计算每个子区域不同粒子间的作用，进而由蛙跳算法继续更新每个粒子的速度；
步骤5)判断是否达到模拟结束条件，判断为是，转至步骤6)，否则转至步骤4)；
步骤6)输出模拟结果；
所述根据更新后的速度与位置重新计算每个子区域不同粒子间的作用，进而由蛙跳算法继续更新每个粒子的速度；具体包括：
针对每个子区域，读取该子区域以及该子区域相邻子区域内的粒子位置信息，采用shift模式，进行不同子区域间的粒子迁移和边界粒子传递；
读取该子区域的一个粒子的位置信息r(t+Δt)，根据更新后的映射关系确定该粒子所在网格，根据所在网格搜寻邻居网格内的其他粒子；
计算该粒子与其他粒子间距离，当距离在截断距离内则构成一个待计算粒子对；
遍历该粒子的每个待计算粒子对，根据粒子对的粒子类型选择相应的粒子作用方式和作用参数，计算对应的作用力或碰撞后的速度更新值；
由该粒子的每个粒子对的作用得到该粒子在(t+Δt)时刻的合力f(t+Δt)；
根据下式得到该粒子在(t+Δt)时刻的速度v(t+Δt)：
所述遍历该粒子的每个待计算粒子对，根据粒子对的粒子类型选择相应的粒子作用方式和作用参数，计算对应的作用力或碰撞后的速度更新值；具体包括：
对于气体粒子间采用拟颗粒作用方式，两个粒子即时更新速度信息；
对于液体粒子间采用光滑粒子动力学作用方式，获得粒子作用力；
对于气液粒子间采用拟颗粒作用方式，即时更新气体拟颗粒的速度信息，并将碰撞过程发生的液体粒子动量变化量除以持续时间转换为作用力；
所述作用参数包括作用截断距离rc、作用特征长度h、黏度系数μ、表面张力系数σ、参考密度ρ0、参考压强p0以及碰撞回复系数e。
2.根据权利要求1所述的基于离散模拟的气液多相系统计算方法，其特征在于，所述模拟设置信息包括模拟区域信息、气体粒子信息、液体粒子信息、气体定向速度、液体注入位置和液体注入速度。
3.根据权利要求1或2所述的基于离散模拟的气液多相系统计算方法，其特征在于，所述采用蛙跳算法同步迭代更新每个粒子的速度与位置，具体包括：
根据粒子在t时刻的速度v(t)、位置r(t)和作用力f(t)，采用蛙跳方式，将速度初步更新至(t+Δt/2)时刻，得到v(t+Δt/2)，从而将位置更新至(t+Δt)时刻得到r(t+Δt)，具体表示如下：
r(t+Δt)＝r(t)+v(t+Δt/2)Δt
其中，m为粒子质量。
4.根据权利要求3所述的基于离散模拟的气液多相系统计算方法，其特征在于，所述更新每个粒子与网格的映射关系；具体包括：
先删除原网格内的映射关系，再依据粒子更新后的位置建立粒子到网格的新映射关系。
5.一种基于离散模拟的气液多相系统计算系统，基于权利要求1-4之一所述的方法实现，其特征在于，所述系统包括：初始设置模块、子区域划分模块、多个并行的处理模块和输出模块；其中，
所述初始设置模块，用于根据模拟设置信息设定待模拟的气液多相系统初始状态；所述初始状态包括气体粒子数目、液体粒子数目、每个气体粒子的初始位置信息及初始速度，每个液体粒子的初始位置信息、初始速度及初始作用力；
所述子区域划分模块，用于对待模拟的气液多相系统进行空间分解得到多个子区域，并分别配置给对应的处理模块；所述处理模块有多个，为并行处理；
所述处理模块，用于对本子区域按照作用截断距离信息进行划分，得到若干个网格；用于根据初始状态，将粒子归属到不同的网格中，建立本子区域每个粒子与网格的映射关系；还用于采用蛙跳算法同步迭代更新每个粒子的速度与位置，进而更新每个粒子与网格的映射关系，再根据更新后的速度与位置重新计算每个子区域不同粒子间的作用，进而由蛙跳算法继续更新每个粒子的速度；并将每次迭代更新的每个粒子的速度与位置发送至输出模块；
所述输出模块，用于输出模拟结果。
CN202110504761.0A 2021-05-10 2021-05-10 一种基于离散模拟的气液多相系统计算方法及系统 Active CN113378445B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110504761.0A CN113378445B (zh) 2021-05-10 2021-05-10 一种基于离散模拟的气液多相系统计算方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110504761.0A CN113378445B (zh) 2021-05-10 2021-05-10 一种基于离散模拟的气液多相系统计算方法及系统

Publications (2)

Publication Number Publication Date
CN113378445A CN113378445A (zh) 2021-09-10
CN113378445B true CN113378445B (zh) 2024-02-02

Family

ID=77571407

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110504761.0A Active CN113378445B (zh) 2021-05-10 2021-05-10 一种基于离散模拟的气液多相系统计算方法及系统

Country Status (1)

Country Link
CN (1) CN113378445B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101727653A (zh) * 2008-10-31 2010-06-09 中国科学院过程工程研究所 一种基于图形处理器的多组分系统离散模拟计算方法
CN102760116A (zh) * 2011-04-29 2012-10-31 中国科学院过程工程研究所 一种基于硬球模型的并行计算方法
CN104143027A (zh) * 2014-08-01 2014-11-12 北京理工大学 一种基于sph算法的流体热运动仿真系统
CN111680456A (zh) * 2020-04-28 2020-09-18 中国科学院深圳先进技术研究院 一种流体力学模拟的方法、装置及存储介质
CN112446942A (zh) * 2020-11-09 2021-03-05 北京达佳互联信息技术有限公司 特效渲染方法、装置、电子设备及存储介质
CN112613243A (zh) * 2020-12-16 2021-04-06 中国科学院深圳先进技术研究院 一种流体力学模拟的方法、装置及计算机可读存储介质

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101727653A (zh) * 2008-10-31 2010-06-09 中国科学院过程工程研究所 一种基于图形处理器的多组分系统离散模拟计算方法
CN102760116A (zh) * 2011-04-29 2012-10-31 中国科学院过程工程研究所 一种基于硬球模型的并行计算方法
CN104143027A (zh) * 2014-08-01 2014-11-12 北京理工大学 一种基于sph算法的流体热运动仿真系统
CN111680456A (zh) * 2020-04-28 2020-09-18 中国科学院深圳先进技术研究院 一种流体力学模拟的方法、装置及存储介质
CN112446942A (zh) * 2020-11-09 2021-03-05 北京达佳互联信息技术有限公司 特效渲染方法、装置、电子设备及存储介质
CN112613243A (zh) * 2020-12-16 2021-04-06 中国科学院深圳先进技术研究院 一种流体力学模拟的方法、装置及计算机可读存储介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"气固两相流内中空多孔催化剂性能的数值模拟";吴秋莹等;《过程工程学报》;第21卷(第7期);775-783 *
单斌.《材料学的纳米尺度计算模拟 从基本原理到算法实现》.华中科技大学出版社,2015,281. *

Also Published As

Publication number Publication date
CN113378445A (zh) 2021-09-10

Similar Documents

Publication Publication Date Title
Ummenhofer et al. Lagrangian fluid simulation with continuous convolutions
Génevaux et al. Simulating Fluid-Solid Interaction.
US7580821B2 (en) Application programming interface for fluid simulations
WO2017092252A1 (zh) 实时虚拟场景中碰撞体之间碰撞探测的方法、终端和存储介质
US9589383B2 (en) Unified position based solver for visual effects
AU2017227323B2 (en) Particle simulation device, particle simulation method, and particle simulation program
Thürey et al. Interactive free surface fluids with the lattice Boltzmann method
CN118013160B (zh) 一种求解多介质流动问题的守恒型界面处理方法及系统
EP3179390A1 (en) Method and apparatus for modeling movement of air bubble based on fluid particles
Yu et al. A robust Delaunay-AFT based parallel method for the generation of large-scale fully constrained meshes
US9613449B2 (en) Method and apparatus for simulating stiff stacks
CN110008387B (zh) 流场可视化实现方法、装置及电子设备
CN113378445B (zh) 一种基于离散模拟的气液多相系统计算方法及系统
KR102123254B1 (ko) 짤리스 엔트로피 및 레이니 엔트로피를 이용한 유체유동 시뮬레이션 방법
Rivara Lepp-bisection algorithms, applications and mathematical properties
EP1979833A2 (en) Object discretization to particles for computer simulation and analysis
Arata et al. Free-form shape modeling by 3D cellular automata
US20200302672A1 (en) Computerized rendering of animated objects using polynomial particle-in-cell method
CN113673186A (zh) 一种基于stl文件的笛卡尔网格快速生成方法
Torres et al. Convex polygon packing based meshing algorithm for modeling of rock and porous media
US20180101975A1 (en) Animating a virtual object
Yoon A study on terrain-surface modeling and searching algorithms for real-time simulation of off-road vehicles
Liu et al. A novel supermesh method for computing solutions to the multi-material Stefan problem with complex deforming interfaces and microstructure
JPH05274277A (ja) 分子動力学計算装置
Zhao et al. INTERACTIVE VISUALIZATION OF LARGE-SCALE 3D SCATTERED DATA FROM A TSUNAMI SIMULATION.

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant