基于SPH 方法一维颗粒仿真分析

程 丽 张 超， 肖道林， 张赵威

（1、沈阳大学机械工程学院，辽宁 沈阳110044 2、中国科学院沈阳自动化研究所，辽宁 沈阳110016）

由于科技的发展和人类对未知领域探索的需要，力学仿真广泛存在于人类科学研究中。但是传统力学仿真所用网格法具有一定局限性。[3-4]

近年，新兴的光滑无网格的方法经过已经成熟，光滑粒子流体力学(SPH)是由核近似和粒子近似两步组成的无网格拉格朗日粒子方法[5]。由于SPH 是在不同时刻计算该区域任意分布粒子上承载的物理量，所以具有极强的自适应性，并且各个物理量均为矢量。基于SPH 方法此种特性，对其一维进行研究可以为高维的工程问题提供参考。

1 SPH 模型

1.1 SPH 基本形式

在SPH 中，通过光滑核函数引进积分表示式[1]。W(xi-xj,h)称为光滑核函数，该函数具有狄拉克函数的性质。其中h 表示光滑长度，? 为积分域，在SPH 中也叫支持域。A（x）j表示物理量A 在坐标xj处的值。x 是空间坐标向量，i 与j 表示粒子序号。

1.2 边界条件的加入

实际工程中流体粒子运动会遇到无法穿透的边界，模拟时为防止粒子穿透边界，采用一种Lennard-Jones 边界力模型

式中，i 和j 分别指边界粒子和内部流体粒子，n1n2为常数，通常取12 和4。D 为最大速度的平方，r0为截止距离，通常取初始距离一半，粒子运动到边界粒子r0范围内将受到边界力的作用。

1.3 时间积分的方法

为了求得粒子在每一个时间步运动状态，引入蛙跳法。根据蛙跳积分的方案，在时间步n 有

2 对此次一维算例进行仿真

图1 0.005s 粒子运动分布状态

图2 0.01975s 粒子运动分布状态

图3 0.025s 粒子分布状态

图4 速度随时间变化曲线

仿真采用颗粒密度1.5kg/m，将每一个粒子看为半径为0.01m。每一个时间步取得到广泛应用的5*10-5s，共计算500 步。仿真结果如下，粒子呈现三角形为向右移动，呈现圆形为向左移动，正方形为边界。

经过1、图2、图3 可以观测到，0.005s 时刻所有粒子还是保持向右移动，0.01975s 开始出现由于边界力作用而向左移动的粒子，最终500 步时即为0.025s 接近一半粒子向左移动，另一半粒子还未运动到边界继续向右移动，图4 对以上加以验证。

3 结论

最终得出SPH 可以模拟一维颗粒状样本变化，面对高维的工程问题只需提升维度如加上环境所需要的外力项与粘性项即可。