[发明专利]一种基于变光滑长度的SPH流体模拟方法

权利要求书

1.一种基于变光滑长度的SPH流体模拟方法，其特征在于包括如下步骤：

步骤一，从对称核近似角度出发修正SPH密度演化方程、动量方程和能量方程，建立一组修正后的变光滑长度SPH方程组；

步骤二，由建立的SPH方程组计算出粒子所受的压力、重力、黏性力，再根据牛顿第二定律计算出合力得到加速度；

步骤三，更新下一步粒子的速度和位置；

步骤四，基于Marching Cubes算法实现流体表面提取，利用GPU纹理缓存的位置、速度和颜色数据进行三维可视化渲染；

Marching Cubes流体表面提取算法涉及网格节点处的密度计算、交点坐标和法矢量计算。采用CUDA并行计算网格节点处的密度值和粒子网格节点密度值，计算思想由影响域改为支持域，即当前网格节点的密度值由邻域范围内所有粒子计算所得，而不是计算一个粒子点对邻域范围内产生的影响。

2.根据权利要求1所述的一种基于变光滑长度的SPH流体模拟方法，其特征在于：所述步骤一具体修正过程如下：

(1)密度方程

传统SPH法中，某点r_i处密度ρ(r_i)的核函数近似式为

其中，W(r_i-r_j,h)称为核函数，h为光滑长度描述核近似的程度，N为i粒子邻近粒子总数。

光滑长度h的选取，应尽量使得粒子的邻近数目在计算中保持不变，以获得全场一致的核近似精度，通常设定粒子i的光滑长度h_i半径区域内包含恒定质量M_sph，即：

其中，d为维数，σ对应一维到三维分别取2、π、公式(2)的引入使得光滑长度在空间分布不一致，并随同粒子的运动变化，即光滑长度同为空间和时间的函数h_i＝h(r_i,t)。

对于变光滑长度核函数近似，采用对称核近似的方法可以得到更高精度的近似值，本发明采用对称光滑长度的方法来对对称函数近似，此时密度核函数近似式为：

其中，

公式(3)两边分别对时间求导数得到：

公式(4)即为用对称核近似的变光滑长度SPH密度演化方程，其中ν_ij＝ν_i-ν_j，在恒定光滑长度时该方程退化为传统SPH密度演化方程。

(2)动量方程

从全守恒SPH方程组思想出发，改用对称的核函数近似来推导动量方程。无耗散流体动力学方程的拉格朗日函数描述如下：

其中，ρ为密度，ν为速度，u为单位质量内能，s为熵。对绝热等熵过程，由热力学第一定律可知：

其中，p为压强。

公式(5)的SPH离散形式为

其中，q＝(r₁,...,r_N,h₁,...,h_N),包含坐标和光滑长度，这一近似过程可理解为把无限维连续介质相空间离散到2N维的有限粒子相空间。

对q的后N个坐标分量—光滑长度引入公式(2)，使所有粒子的光滑区域内包含恒定质量物质。其约束条件描述如下：

采用带约束的拉格朗日方程：

令由公式(9)前N个方程以及梯度核函数可得

公式(10)即为采用对称核近似的变光滑长度SPH动量方程，f_i为修正系数。

(3)能量方程

由热力学第一定律可知

其中，T为温度，s为熵，u为单位质量内能，p为压强，v为比体积。若没有热源产生，方程可以简化为：

公式(12)的SPH离散形式为采用对称核近似的变光滑长度SPH能量方程。

。

3.根据权利要求1所述的一种基于变光滑长度的SPH流体模拟方法，其特征在于：所述步骤二为有效计算压力，将不可压缩流体视为若弱可压缩流体。引入弱可压状态方程：

其中，参数P₀为参考压强，γ为常数。P₀和γ共同用于控制计算中流体密度，以保证模拟流场的不可压缩性。