CN104537175A - 一种基于sph算法的流体模拟方法及装置 - Google Patents

Abstract

本发明适用于计算机仿真技术领域，提供了一种基于SPH算法的流体模拟方法及装置，包括：初始化流场，生成离散的粒子；计算当前流场中每个所述粒子的密度；根据所述粒子的密度计算每个所述粒子的物理属性，所述物理属性包括所述粒子的涡流约束；根据每个所述粒子的物理属性计算每个所述粒子的速度变化率；根据每个所述粒子的速度变化率更新每个所述粒子的速度和位置。本发明通过将涡流约束无缝结合到SPH算法中，增强了基于SPH算法模拟出的流体的细节表现，提高了流体模拟的保真度。

Description

一种基于SPH算法的流体模拟方法及装置

技术领域

本发明属于计算机仿真技术领域，尤其涉及一种基于SPH算法的流体模拟方法及装置。

背景技术

在计算机图形学领域，近年来，基于粒子的拉格朗日方法逐渐成为实现流体模拟的主要工具，特别是平滑粒子流体动力学(Smoothed PartideHydrodynamics，SPH)方法，其能够更为自然地描述介质运动的历史过程，模拟更多的流体细节，例如泡沫、水花等，以及处理更为复杂的流体表面，且该方法计算量小，易于实现。

然而，由于粒子的离散化以及模拟区域的限制，SPH方法不可避免地增加了数值耗散，这将导致流体模拟过程中很多微小细节的损失。

发明内容

本发明实施例的目的在于提供一种基于SPH算法的流体模拟方法及装置，旨在解决现有技术中基于SPH算法的流体模拟方法导致很多微小细节损失的问题。

本发明实施例是这样实现的，一种基于SPH算法的流体模拟方法，包括：

初始化流场，生成离散的粒子；

计算当前流场中每个所述粒子的密度；

根据所述粒子的密度计算每个所述粒子的物理属性，所述物理属性包括所述粒子的涡流约束；

根据每个所述粒子的物理属性计算每个所述粒子的速度变化率；

根据每个所述粒子的速度变化率更新每个所述粒子的速度和位置。

本发明实施例的另一目的在于提供一种基于SPH算法的流体模拟装置，包括：

初始化单元，用于初始化流场，生成离散的粒子；

密度计算单元，用于计算当前流场中每个所述粒子的密度；

物理属性计算单元，用于根据所述粒子的密度计算每个所述粒子的物理属性，所述物理属性包括所述粒子的涡流约束；

速度变化率计算单元，用于根据每个所述粒子的物理属性计算每个所述粒子的速度变化率；

模拟单元，用于根据每个所述粒子的速度变化率更新每个所述粒子的速度和位置。

本发明实施例通过将涡流约束无缝结合到SPH算法中，增强了基于SPH算法模拟出的流体的细节表现，提高了流体模拟的保真度。

附图说明

图1是本发明实施例提供的基于SPH算法的流体模拟方法的实现流程图；

图2是本发明实施例提供的SPH算法原理示意图；

图3是本发明实施例提供的基于SPH算法的流体模拟方法S103的具体实现流程图；

图4是本发明实施例提供的基于SPH算法的流体模拟方法所模拟的波浪效果对比图；

图5是本发明实施例提供的基于SPH算法的流体模拟方法所模拟的流体与固体交互的场景效果对比图；

图6是本发明实施例提供的基于SPH算法的流体模拟装置的结构框图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

对流体进行模拟，无论是采用基于网格的欧拉法还是采用基于粒子的拉格朗日法，多数都是以Navier-stoke方程作为基本的流体运动物理方程进行力的求解与模拟。所述Navier-Stokes方程，是描述不可压缩流体的动量守恒运动方程，对于恒温、带有滞粘性的流体，其基本的物理量包括速度、质量以及压强。

Navier-Stokes方程的通常表示形式为：

$\frac{\∂\stackrel{\→}{u}}{\∂t}+\stackrel{\→}{u}\·\▿\stackrel{\→}{u}+\frac{1}{\mathrm{\ρ}}\▿p=\stackrel{\→}{g}+v\▿\·\▿\stackrel{\→}{u}---\left(1\right)$

$\▿\·\stackrel{\→}{u}=0---\left(2\right)$

其中，所述表示速度矢量；所述ρ表示粒子的密度；所述p表示粒子的压强；所述代表重力加速度，其值在三维空间中约为(0，0，-9.8)m/s³，且可以不局限于表示重力加速度，任何作用在整个流体之上的作用力都可以用这一项来表示；所述v是运动黏度，其可以用来衡量一种流体有多黏，例如，糖浆有非常高的黏度，而酒精的黏度就很低，运动黏度可以用来测量流体流动时形变所产生的阻力。

在本发明实施例中，采用SPH方法来对流体进行模拟，由于粒子的数量是固定的并且每个粒子的质量是固定的，质量守恒完全可以得到保证，因此，式(2)可以完全忽略，而式(1)左侧表达式的前两项可以用随体导数代替。随着粒子运动代替了流体运动，无需再构建网格来求解偏微分方程，速度场的随体导数简化了粒子速度对时间的导数，这意味着对流项在粒子系统中已经不再需要。因此Navier-Stokes方程可简化为：

$\mathrm{\ρ}\frac{D\stackrel{\→}{u}}{\mathrm{Dt}}=\mathrm{\ρ}\stackrel{\→}{g}+v\▿\·\▿\stackrel{\→}{u}-\▿p---\left(3\right)$

其中，式(3)右边的每一项分别对应重力，粘性力及流体内部压力，所述重力即为常数-9.8m/s³。SPH算法即通过求解式(3)来得到每个粒子的总加速度，最终完成对粒子位置的更新。

图1示出了本发明实施例提供的基于SPH算法的流体模拟方法的实现流程，详述如下：

在S101中，初始化流场，生成离散的粒子。

SPH方法通过离散的采样点来计算连续的流场，其把整个被模拟的流体离散成采样点，这些离散的采样点被称为光滑粒子。在SPH方法中，粒子可以分布于被模拟空间的任意位置，每个粒子都被认为占据着流体空间的一个片段，且为了得到更准确的加权量，采样的粒子必须比较密集。

在本实施例中，可以通过C++或者OpengL来实现流场的初始化。

在S102中，计算当前流场中每个所述粒子的密度。

SPH方法本质上是一种插值方法，其基本思想是通过插值周围粒子的属性来计算空间中任意位置r处具有连续性的属性A，即：

A(r)＝∫A(r′)W(r-r′,h)dr′   (4)

其中，W()是以h为核半径的平滑核函数，所述核半径h是一个可调的参数，用于改变流体的光滑程度。

将(4)式的积分近似地离散化成数值的求和形式：

$A\left(r\right)={\mathrm{\Σ}}_j{A}_j\frac{m_j}{{\mathrm{\ρ}}_j}W(r-r_j,h)---\left(5\right)$

其中，所述j为所述当前粒子i的邻域粒子j，所述m_j为所述粒子j的质量，所述ρ_j是邻域粒子j的密度，所述r_j为邻域粒子j的位置，A_j为邻域粒子j在其位置r_j处的属性。

图2示出了本发明实施例提供的SPH算法原理示意图。

则根据式(5)，当A_j属性为密度时，即A_j＝ρ_j时，得到当前粒子i的密度计算公式为：

${\mathrm{\ρ}}_j={\mathrm{\Σ}}_j{\mathrm{\ρ}}_j\frac{m_j}{{\mathrm{\ρ}}_j}W(r-r_j,h)={\mathrm{\Σ}}_j{m}_{j}W(r_i-r_j,h)---\left(6\right)$

在S103中，根据所述粒子的密度计算每个所述粒子的物理属性，所述物理属性包括所述粒子的涡流约束。

每个粒子都包含一些物理属性，其中有三个基本物理属性：粒子的质量、粒子的半径和粒子的位置，除了以上三个基本属性之外，还可以根据模拟需求附加上其他的物理属性，所述物理属性包括粒子的压力、粘性力、重力，等等。

在S102计算得到每个粒子的密度后，便可以进一步计算得到每个粒子的其余物理属性，其中：

当前粒子i的压力的计算公式如下：

$F_i^{\mathrm{pressure}}=-\▿p\left(r_i\right)=-{\mathrm{\Σ}}_j{m}_j\left(\frac{p_i+p_j}{2{\mathrm{\ρ}}_j}\right)\▿W(r_i-r_j,h)---\left(7\right)$

当前粒子i的粘性力的计算公式如下：

$F_i^{\mathrm{vis}\cos \mathrm{ity}}=v\▿\·\▿\stackrel{\→}{u}\left(r_i\right)=v{\mathrm{\Σ}}_j{m}_j\left(\frac{{\stackrel{\→}{u}}_j-{\stackrel{\→}{u}}_i}{{\mathrm{\ρ}}_j}\right){\▿}^{2}W(r_i-r_j,h)---\left(8\right)$

由于SPH算法本身存在的非物理数值耗散减弱了一些微小的流动特征，因此，在本实施例中，为了能够更加逼真地模拟流体，除了模拟粒子的压力和粘性力，还需要将涡流约束结合到SPH算法中。所述涡流约束，最早用于对直升机周围复杂混乱流动域的数值进行计算，之后被用来模拟烟雾的湍流现象，然而，目前的涡流约束都是基于网格法来进行模拟的，在本实施例中，需要基于粒子法来实现对涡流约束的模拟。

具体地，如图3所示：

在S301中，计算每个所述粒子的涡量。

将涡量公式转换成SPH格式，得到：

${\mathrm{\ω}}_i=\▿\⊗\stackrel{\→}{u}={\mathrm{\Σ}}_j{m}_j{\stackrel{\→}{u}}_{\mathrm{ij}}\⊗\▿W(r_i-r_j,h)---\left(9\right)$

其中，所述所述为所述当前粒子i的速度，所述为所述邻域粒子j的速度，所述W()是以h为核半径的平滑核函数，所述h为一常数，所述r_i为所述当前粒子i的位置，所述r_j为所述邻域粒子j的位置，是叉乘运算符。

式(9)计算出的涡量提供小规模结构，每一片小涡量可以被认为是一个朝特定方向旋转的流场。

在S302中，根据每个所述粒子的涡量，为每个所述粒子构建校正力场。

在得到每一个粒子的涡量后，便可以通过下式为每一个粒子构建一个校正力场来防止小尺度细节的消失：

$F_i^{\mathrm{vorticity}}=\∈(N\⊗{\mathrm{\ω}}_i)---\left(10\right)$

其中，所述为当前粒子i的校正力场，所述∈为一常数，用于控制涡流约束的尺度，所述N为指向涡流旋转中心的归一化向量，N＝η/|η|，而所述将式(10)转换成SPH格式，得到：

${\mathrm{\η}}_i={\mathrm{\Σ}}_j(m_j/{\mathrm{\ρ}}_j)\left|{\mathrm{\ω}}_j\right|\▿W(r_i-r_j,h)---\left(11\right)$

在S104中，根据每个所述粒子的物理属性计算每个所述粒子的速度变化率。

在S105中，根据每个所述粒子的速度变化率更新每个所述粒子的速度和位置。

在确定了粒子的物理属性之后，根据粒子的物理属性对粒子加速度的大小和方向的影响，基于SPH算法计算每个粒子的速度变化率，并进一步确定每个粒子下一个时间单位的速度和位置，如此循环执行S102至S105，即完成了对流体的模拟。

图4示出了基于上文实施例所述的基于SPH算法的流体模拟方法所模拟的波浪效果对比图，在实验中，为了模拟波浪击打的效果，在模拟的区域空间中设计一个周期性摆的虚拟平面板。在图4中，上图为传统SPH算法模拟的结果，下图是基于本发明的模拟结果，且两种模拟方法均采用了相同的参数设置，从图4中可以看出，下图相对于上图在流体表面上产生了小尺度的波纹，还原了流体表面的微小细节，模拟的保真度高。

图5示出了基于上文实施例所述的基于SPH算法的流体模拟方法所模拟的流体与固体交互的场景效果对比图，在实验中，在流体空间内放入一个周期性运动的固体，并实现了固体与流体粒子的碰撞检测运算。在图5中，上图为传统SPH算法模拟的结果，下图是基于本发明的模拟结果，且两种模拟方法均采用了相同的参数设置，从图5中可以看出，固体在流体空间运动的过程中，传统SPH算法模拟的流体表面比较平坦，而本发明模拟的流体表面的扰动效果明显，模拟效果细节更为丰富。

本发明实施例通过将涡流约束无缝结合到SPH算法中，增强了基于SPH算法模拟出的流体的细节表现，提高了流体模拟的保真度。

对应于上文实施例所述的基于SPH算法的流体模拟方法，图6示出了本发明实施例提供的基于SPH算法的流体模拟装置的结构框图，为了便于说明，仅示出了与本实施例相关的部分。

参照图6，该装置包括：

初始化单元61，初始化流场，生成离散的粒子。

密度计算单元62，计算当前流场中每个所述粒子的密度。

物理属性计算单元63，根据所述粒子的密度计算每个所述粒子的物理属性，所述物理属性包括所述粒子的涡流约束。

速度变化率计算单元64，根据每个所述粒子的物理属性计算每个所述粒子的速度变化率。

模拟单元65，根据每个所述粒子的速度变化率更新每个所述粒子的速度和位置。

可选地，所述物理属性计算单元63包括：

涡量计算子单元，计算每个所述粒子的涡量。

构建子单元，根据每个所述粒子的涡量，为每个所述粒子构建校正力场，以得到所述粒子的涡流约束。

可选地，所述涡量计算子单元具体用于：

根据计算每个所述粒子的涡量，其中，所述ω_i为当前粒子i的涡量，所述j为所述当前粒子i的邻域粒子j，所述m_j为所述粒子j的质量，所述所述为所述当前粒子i的速度，所述为所述邻域粒子j的速度，所述W()是以h为核半径的平滑核函数，所述h为一常数，所述r_i为所述当前粒子i的位置，所述r_j为所述邻域粒子j的位置；

所述构建子单元具体用于：

根据构建每个所述粒子的校正力场，其中，所述为所述当前粒子i的校正力场，所述∈为一常数，所述N为指向涡流旋转中心的归一化向量，N＝η/|η|，所述

将所述 $F_i^{\mathrm{vorticity}}=\∈(N\×\⊗{\mathrm{\ω}}_i)$ 转换成 ${\mathrm{\η}}_i={\mathrm{\Σ}}_j(m_j/{\mathrm{\ρ}}_j)\left|{\mathrm{\ω}}_j\right|\▿W(r_i-r_j,h),$ 其中，所述ρ_j为所述邻域粒子j的密度，所述ω_j为所述邻域粒子j的涡量。

可选地，所述物理属性还包括所述粒子的压力，所述物理属性计算单元63具体用于：

根据 $F_i^{\mathrm{pressure}}=-{\mathrm{\Σ}}_j{m}_j\left(\frac{p_i+p_j}{2{\mathrm{\ρ}}_j}\right)\▿W(r_i-r_j,h)$ 计算所述当前粒子i的压力其中，所述p_i为所述当前粒子i的压强，所述p_j为所述邻域粒子j的压强。

可选地，所述物理属性还包括所述粒子的粘性力，所述物理属性计算单元63具体用于：

根据 $F_i^{\mathrm{vis}\cos \mathrm{ity}}=v{\mathrm{\Σ}}_j{m}_j\left(\frac{{\stackrel{\→}{u}}_j-{\stackrel{\→}{u}}_i}{{\mathrm{\ρ}}_j}\right){\▿}^{2}W(r_i-r_j,h)$ 计算所述当前粒子i的粘性力其中，所述v是运动黏度。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于平滑粒子流体动力学SPH算法的流体模拟方法，其特征在于，包括：

初始化流场，生成离散的粒子；

计算当前流场中每个所述粒子的密度；

根据所述粒子的密度计算每个所述粒子的物理属性，所述物理属性包括所述粒子的涡流约束；

根据每个所述粒子的物理属性计算每个所述粒子的速度变化率；

根据每个所述粒子的速度变化率更新每个所述粒子的速度和位置。

2.如权利要求1所述的方法，其特征在于，所述根据所述粒子的密度计算每个所述粒子的物理属性包括：

计算每个所述粒子的涡量；

根据每个所述粒子的涡量，为每个所述粒子构建校正力场，以得到所述粒子的涡流约束。

3.如权利要求2所述的方法，其特征在于，所述计算每个所述粒子的涡量包括：

根据计算每个所述粒子的涡量，其中，所述ω_i为当前粒子i的涡量，所述j是为所述当前粒子i的邻域粒子j个数，所述m_j为所述邻域粒子j的质量，所述所述为所述当前粒子i的速度，所述为所述邻域粒子j的速度，所述W()是以h为核半径的平滑核函数，所述h为一常数，所述r_i为所述当前粒子i的位置，所述r_j为所述邻域粒子j的位置；

所述根据每个所述粒子的涡量，为每个所述粒子构建校正力场，以得到所述粒子的涡流约束包括：

根据构建每个所述粒子的校正力场，其中，所述为所述当前粒子i的校正力场，所述∈为一常数，所述N为指向涡流旋转中心的归一化向量，N＝η/|η|，所述

将所述 $F_i^{\mathrm{vorticity}}=\∈(N\⊗{\mathrm{\ω}}_i)$ 转换成 ${\mathrm{\η}}_i={\mathrm{\Σ}}_j(m_j/{\mathrm{\ρ}}_j)\left|{\mathrm{\ω}}_j\right|\▿W(r_i-r_j,h),$ 其中，所述ρ_j为所述邻域粒子j的密度，所述ω_j为所述邻域粒子j的涡量。

4.如权利要求3所述的方法，其特征在于，所述物理属性还包括所述粒子的压力，所述根据所述粒子的密度计算每个所述粒子的物理属性包括：

根据 $F_i^{\mathrm{pressure}}=-{\mathrm{\Σ}}_j{m}_j\left(\frac{p_i+p_j}{2{\mathrm{\ρ}}_j}\right)\▿W(r_i-r_j,h)$ 计算所述当前粒子i的压力其中，所述p_i为所述当前粒子i的压强，所述p_j为所述邻域粒子j的压强。

5.如权利要求3所述的方法，其特征在于，所述物理属性还包括所述粒子的粘性力，所述根据所述粒子的密度计算每个所述粒子的物理属性包括：

根据 $F_i^{\mathrm{vis}\cos \mathrm{ity}}=v{\mathrm{\Σ}}_j{m}_j\left(\frac{{\stackrel{\→}{u}}_j-{\stackrel{\→}{u}}_i}{{\mathrm{\ρ}}_j}\right){\▿}^{2}W(r_i-r_j,h)$ 计算所述当前粒子i的粘性力所述v为运动黏度。

6.一种基于平滑粒子流体动力学SPH算法的流体模拟装置，其特征在于，包括：

初始化单元，用于初始化流场，生成离散的粒子；

密度计算单元，用于计算当前流场中每个所述粒子的密度；

物理属性计算单元，用于根据所述粒子的密度计算每个所述粒子的物理属性，所述物理属性包括所述粒子的涡流约束；

速度变化率计算单元，用于根据每个所述粒子的物理属性计算每个所述粒子的速度变化率；

模拟单元，用于根据每个所述粒子的速度变化率更新每个所述粒子的速度和位置。

7.如权利要求6所述的装置，其特征在于，所述物理属性计算单元包括：

涡量计算子单元，用于计算每个所述粒子的涡量；

构建子单元，用于根据每个所述粒子的涡量，为每个所述粒子构建校正力场，以得到所述粒子的涡流约束。

8.如权利要求7所述的装置，其特征在于，所述涡量计算子单元具体用于：

根据计算每个所述粒子的涡量，其中，所述ω_i为当前粒子i的涡量，所述j为所述当前粒子i的邻域粒子j，所述m_j为所述粒子j的质量，所述所述为所述当前粒子i的速度，所述为所述邻域粒子j的速度，所述W()是以h为核半径的平滑核函数，所述h为一常数，所述r_i为所述当前粒子i的位置，所述r_j为所述邻域粒子j的位置；

所述构建子单元具体用于：

根据构建每个所述粒子的校正力场，其中，所述为所述当前粒子i的校正力场，所述∈为一常数，所述N为指向涡流旋转中心的归一化向量，N＝η/|η|，所述

将所述 $F_i^{\mathrm{vorticity}}=\∈(N\⊗{\mathrm{\ω}}_i)$ 转换成 ${\mathrm{\η}}_i={\mathrm{\Σ}}_j(m_j/{\mathrm{\ρ}}_j)\left|{\mathrm{\ω}}_j\right|\▿W(r_i-r_j,h),$ 其中，所述ρ_j为所述邻域粒子j的密度，所述ω_j为所述邻域粒子j的涡量。

9.如权利要求8所述的装置，其特征在于，所述物理属性还包括所述粒子的压力，所述物理属性计算单元具体用于：

根据 $F_i^{\mathrm{pressure}}=-{\mathrm{\Σ}}_j{m}_j\left(\frac{p_i+p_j}{2{\mathrm{\ρ}}_j}\right)\▿W(r_i-r_j,h)$ 计算所述当前粒子i的压力其中，所述p_i为所述当前粒子i的压强，所述p_j为所述邻域粒子j的压强。

10.如权利要求8所述的装置，其特征在于，所述物理属性还包括所述粒子的粘性力，所述物理属性计算单元具体用于：

根据 $F_i^{\mathrm{vis}\cos \mathrm{ity}}=v{\mathrm{\Σ}}_j{m}_j\left(\frac{{\stackrel{\→}{u}}_j-{\stackrel{\→}{u}}_i}{{\mathrm{\ρ}}_j}\right){\▿}^{2}W(r_i-r_j,h)$ 计算所述当前粒子i的粘性力其中，所述v是运动黏度。