一种基于密度分离式SPH模型的多相流模拟方法

技术领域

本发明属于计算流体力学领域，具体涉及一种基于密度分离式SPH(SmoothedParticle Hydrodynamics,光滑粒子流体动力学)模型的多相流模拟方法。

背景技术

传统的数值模拟方法大多使用欧拉型网格方法，如有限差分法(FDM)、有限体积法(FVM)和有限元法(FEM)，来近似基于不规则网格插值的Navier-Stokes动量方程。然而，传统的网格方法需要特殊的方法来处理自由表面流动问题。例如，有限元法(FEM)不能有效地解决大的流体变形问题，需要调整网格，有限体积法(FDM)需要跟踪自由表面的变化，即使使用拉格朗日动态网格，也可能出现网格扭曲的问题。

相比网格法，无网格法在模拟具有自由表面的不可压缩流、边界变形和移动界面中具有特殊优势。作为无网格方法的典型代表，SPH的应用领域广泛，包括弹性流、多相流、准不可压流、重力流、多孔介质流、热传导、冲击模拟、热传递和质量流动等问题。

然而，当模拟涉及跨越介质的传播与交换的问题时，由于密度、速度等物理场在相界面处的不连续，SPH多相模型比单相模型更具复杂性和挑战性。为此，学者们针对多相界面处理提出并发展了许多方法，如密度重整、表面张力模型、界面力模型、体积求和法等，但现有的多相流SPH模型仍然存在很多不足，例如精度低、稳定性差、计算时间长，亟需一种能够同时兼顾数值稳定性和数值精度的多相流SPH模型。

发明内容

针对现有技术中存在的“SPH多相流模型精度不足、计算时间长”的问题，本申请提出一种基于密度分离式SPH模型的多相流模拟方法，通过在求解连续性方程时分离不同相的粒子，同时在求解动量方程时，使用修正的压力值代替异相粒子的压力，达到提升计算精度和稳定性，允许在多相流计算时采用更大的CFL数，进而提升计算效率的目的。

为解决上述技术问题，本发明采取的技术方案是：

一种基于密度分离式SPH模型的多相流模拟方法，包括以下步骤：

1)依据待模拟的多相流动过程开始时刻各相流体的分布情况及具体尺寸，用不同类型的SPH粒子离散不同相的流体，建立相应的多相流数值模型；在某一时间步内，遍历所有粒子，找到所有与粒子i发生相互作用的粒子j，

2)如果粒子i属于流体粒子，那么使用经SPH法离散后的连续性方程计算粒子i的密度增量，只考虑支持域内的虚粒子或与粒子i同相的流体粒子对粒子i密度增量的贡献，

3)如果粒子i属于流体粒子，使用经SPH法离散后的Navier-Stokes方程计算粒子i的速度增量，考虑支持域内所有相的流体粒子和虚粒子对粒子i速度变化的贡献，但是计算与粒子i不同相的流体粒子j对粒子i速度变化的贡献时，使用粒子j所在位置处压力的Shepard插值结果代替粒子j本身的压力参与计算，

4)如果i属于边界粒子，使用Shepard插值求解粒子i的压力和速度；粒子i的压力插值中，支持域内所有粒子都参与插值计算：

上式中，f代表流体粒子，p

同一个边界粒子i应当具有多个不同的速度以参与不同相的流体粒子的计算，

其中，u

例如，当模拟气-液两相流问题时，边界粒子i应当具有两个不同的速度u

5)对所有流体粒子，按步骤2)和3)求得粒子的速度增量和密度增量后，依照计算所用的时间推进格式，更新粒子的密度和速度，进入下一时间步的求解。

根据本发明的优选方案，步骤2)中，判断粒子i和j发生相互作用的依据为：粒子i和j的距离r

根据本发明的优选方案，所述步骤2)中，具体按下式计算粒子i的密度增量：

其中,ρ

根据本发明的优选方案，在步骤2)中，可以直接离散连续性方程求解粒子i的密度增量，也可在连续性方程中添加密度耗散项求解，具体为：在连续性方程中加入密度耗散项后的连续性方程如下：

式中，

根据本发明的优选方案，所述步骤3)中，具体按下式计算粒子i的速度增量：

其中，

π

上式中，Ω

根据本发明的优选方案，步骤4)中，可以在边界粒子所在位置处进行Shepard插值，也可使用映射点方法，从边界向内沿边界的法向延伸半个粒子间距作为映射点，使用映射点处压力和速度的Shepard插值结果作为边界粒子的速度和压力值。

对任意粒子i，无论其是流体粒子或边界粒子，压力p

本申请的有益效果为：在多相流相界面处求解不同相粒子相互作用时，只考虑同相粒子对某一粒子密度的影响，认为异相粒子对密度没有影响，使数值模型更符合物理实际，提升模型精度；考虑异相粒子对速度的影响时，以异相粒子所在位置压力的Shepard插值代替异相粒子的压力，从而提升相界面处的数值稳定性。现有技术中的SPH多相流格式由于精度低、数值稳定性差，往往需要极小的CFL数(0.05-0.1)和时间步长以保证数值稳定，本申请所提出的新型分离式SPH多相流模型能够在保证数值精度的情况下，以较大的CFL条件数(0.2)计算仍然能够确保数值计算的稳定性，进而大幅提升计算效率，节约计算成本。

附图说明

图1对处于相界面附近的某一个中心粒子，(a)所有与其同相粒子形成的集合χ示意图，(b)支持域内所有粒子形成的集合Ω示意图；

图2相界面附近边界粒子压力和速度插值示意图，(a)压力插值示意图，(b)速度插值示意图；

图3考虑气-液两相的数值波浪水槽模型布置示意图；

图4模拟中粒子布置图；

图5相对波高H

图6相对波高H

图7孤立波破碎模拟示意图。

图8研究不同坡度时的破碎指标和破碎类型时，本方法预测结果与经验公式结果对比。

图9相对波高H

具体实施方式

下面结合具体实施方式对本发明做进一步阐述和说明。所述实施例仅是本公开内容的示范且不圈定限制范围。本发明中各个实施方式的技术特征在没有相互冲突的前提下，均可进行相应组合。

实例1、下面以考虑气-液两相流的数值波浪水槽造波为例，来说明本发明的技术方案。模拟包括以下步骤：

根据如附图3所示的计算模型，建立数值水槽，同时初始化模拟所需要的各项参数。数值模拟中监测波浪高度的位置距离造波板x＝0.4m，水深d＝0.114m，相对波高H

1)开始计算后，每一步的计算步骤为：

1a)在某一时间步内，遍历所有粒子，找到所有与粒子i发生相互作用的粒子j，判断粒子i和j发生相互作用的依据为：粒子i和j的距离r

1b)如果i属于流体粒子，使用经SPH法离散后的连续性方程计算粒子i的密度增量，但只考虑支持域内的虚粒子或与粒子i同相的流体粒子对粒子i密度增量的贡献；使用经SPH法离散后的Navier-Stokes方程计算粒子i的速度增量，考虑支持域内所有相的流体粒子和虚粒子对粒子i速度变化的贡献，但是计算与粒子i不同相的流体粒子j对粒子i速度变化的贡献时，使用粒子j所在位置处压力的Shepard插值结果代替粒子j本身的压力参与计算；

1c)如果i属于边界粒子，使用Shepard插值求解粒子i的压力和速度。1d)对所有流体粒子，按步骤2b)求得粒子的速度增量和密度增量后，依照计算所用的时间推进格式(如预测-校正时间格式)，更新粒子的密度和速度，进入下一时间步的求解。

2)达到计算所需要的物理时间后，停止计算，输出结果。

附图5给出了相对波高H

实例2、当一个孤立波在近岸传播并遇到斜坡时，波浪在爬坡过程中会由于当地水深的大小变化，逐渐发生浅水变形和破碎。在破碎过程中，由于伴随着相界面的剧烈变形和气-液两相之间剧烈的掺混，现有的多相流数值模拟方法往往难以准确模拟，或者需要极小的时间步长以保证模拟过程中的数值稳定性。下面以考虑气-液两相流的孤立波破碎过程为例，来说明本发明的技术方案。模拟包括以下步骤：

根据如附图7所示的计算模型，建立数值水槽，同时初始化模拟所需要的各项参数。数值模拟中造波板距离斜坡坡脚1.0m，水深d＝0.2m，相对波高H

1)开始计算后，每一步的计算步骤为：

1a)在某一时间步内，遍历所有粒子，找到所有与粒子i发生相互作用的粒子j，判断粒子i和j发生相互作用的依据为：粒子i和j的距离r

1b)如果i属于流体粒子，使用经SPH法离散后的连续性方程计算粒子i的密度增量，但只考虑支持域内的虚粒子或与粒子i同相的流体粒子对粒子i密度增量的贡献；使用经SPH法离散后的Navier-Stokes方程计算粒子i的速度增量，考虑支持域内所有相的流体粒子和虚粒子对粒子i速度变化的贡献，但是计算与粒子i不同相的流体粒子j对粒子i速度变化的贡献时，使用粒子j所在位置处压力的Shepard插值结果代替粒子j本身的压力参与计算；

1c)如果i属于边界粒子，使用Shepard插值求解粒子i的压力和速度。1d)对所有流体粒子，按步骤2b)求得粒子的速度增量和密度增量后，依照计算所用的时间推进格式(如预测-校正时间格式)，更新粒子的密度和速度，进入下一时间步的求解。

2)达到计算所需要的物理时间后，停止计算，输出结果。

通常而言，如果斜坡越缓或孤立波波高越高，也就是当斜坡坡度s减少和相对波高H

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。