３数值试验及结果分析
３．１ 参数设置及粒子初始化 数值试验表明，粒子尺寸较大时，计算结果易产生
振荡，而粒子尺寸太小，粒子数量增多，计算耗时过长． 根据实际计算条件选为ｄｏ＝２．５ ｍｒｆｌ．２组实验的固边 界粒子均为１ ７２２个；水柱高度为１８ ｃｍ的水粒子为 ５ １８４个，２４ ｃｍ的水粒子６ ９１２个．粒子初始布置如图 ３所示．
粒子被判定为自由表面粒子（见图１）．
叫哔兰
图１边界粒子布置 Ｆｉｇ．１ ＣｏＩＩｆｉｇｕ豫６０ｎ ０ｆ ｒｉｇｍ ｂ０删Ｉｄａｒｙ ｐａｒｔｉｃｌ豁
１．７摩阻项 考虑到水柱倒塌过程中水槽底部与侧壁的摩阻效
应，方程（２）右端添加摩阻力项Ｆｉ，可得
掣：一上ｖｐ＋ｙ铲髓＋Ｆ＋Ｂ
（１９）
Ｑｆ
ｐ
采用阻力平方公式
度满足
ｎ？＜胁。
（１８）
即被判定为自由表面粒子，其中ＪＢ为小于１的数．本研
究大量数值实验表明，卢＝ｏ．９９较适宜．对于自由表面
粒子，其压力直接等于大气压，无需参与压力计算．
１．６固边界处理
固边界由大小与流体粒子相同的边界粒子组成，
它们与流体粒子一起参与压力计算，但不进行速度和
坐标校正；边界粒子以外另设了两层虚粒子，以免边界
扩散模型是用来计算空间二阶导数项的，亦称
Ｌａｐｌａｃｅ模型．由物理量厂的扩散问题
誓＝矿町
的解析解，得到出时间内周围粒子对ｉ的物理量输运
量为
馘ｚ∑｛皤ｉ一馘叫｝＝
半去军奶圳删。＿Ｉ）（１２）
因为式（１１）的解析解是在无限空间中按高斯函
数分布求得的，而核子函数仅限于有限的加权半径范
围内，所以引入Ａ来校正这种误差，即
Ｆｉｇ．５
实测结果对比（第２组次） Ｃｏｍｐａｒｉｓｏｎ ｂｅｔｗｅｅｎ ｓｉｍｕｌａｔｅｄ ａｎｄ ｍｅａｓｕｒｅｄ ｒｅｓｕｌｔｓ ｏｆ ｐｒｅｓｓｕｒｅ ｐｒｏｃｅｓｓ ｏｎ ｒｉｇｈｔ ｂｏａｒｄ
ａｆｔｅｒ ｗａｔｅｒ ｃｏｌｕｍｎ ｃｏｌｌａｐｓｅ（Ｃａｓｅ ２）
万方数据
Ｆｉｇ．６
图６ Ａ、Ｂ时刻流场图（第１组次） Ｆｌｏｗ ｆｉｅｌｄｓ ａｔ ｔｈｅ ｔｉｍｅ ｏｆ ｔｈｅ ｔｗｏ ｐｒｅｓｓｕｒｅ ｐｅａｋｓＡ ａｎｄ Ｂ（Ｃａｓｅ １）
·１０３４·
天津大学学报
第３９卷第９期
衄，·：一垒ＶＰ川
（２４）
Ｐ
比∥＝留“＋Ａｕｉ一
（２５）
∥１＝，？＋丢（∥１＋Ｈ７）ａｔ
（２６）
２水柱倒塌物理试验
该试验是在一个自制有机玻璃水槽中进行的（图 ２）．水槽右侧为挡板，沿挡板中心轴线装有１０个压力
传感器（图２（ｂ）），用于测定水柱倒塌后对挡板的冲 击压力．在一个由活动挡板和水槽壁形成的密封空间 内预装一定量清水．试验开始时，扣动扳机，活动挡板 在重锤的带动下，在极短时间内（小于Ｏ．０５ Ｓ）完成向 上翻转，水柱倒塌过程基本不受挡板影响．
图５（ａ）给出了第２组２４ ｃｍ水柱当黏滞系数取 不同值时挡板压力过程的计算结果．计算中因未加入 摩阻项，所以峰点相位与实测结果相比均有一定提前． 从各组计算结果看，黏滞系数对压力峰的尖锐程度有 显著影响．黏滞系数取矽＝１．０１×ｉ０山Ｉｎ２／ｓ时，峰型 最尖锐．随着系数增大，峰型逐渐钝化．当黏滞系数取 ／／＝３．５０×１０。６ Ｉｎ２／ｓ时，峰型与实测结果最接近．
（即ｉ＝竿∑∽一ｚ）埘（１０一，ｉ ｌ） （１６）
式中
Ａｉ＝∑训（１ ｏ—ｔ １）ｌ，Ｊ—ｔ ２
（１７）
１．５流体自由表面的确定 流体自由表面由位于自由表面的粒子形成．由于
自由表面粒子周围的粒子数量较少，其粒子数密度低 于内部粒子的水平，这样就可以通过粒子数密度的大
万方数据
小来判定粒子是否是自由表面粒子．若瞬时粒子数密
采用北京水科院研制的ＤＪ８００数据采集系统，每 个传感器采集时间间隔为０．００５ Ｓ，每隔０．０５ Ｓ １０个 传感器轮流采集１次．水柱宽度为１８ ｃｍ，高度分为 １８ ｃｍ和２４ ｃｍ ２个组次．
（ａ）实验水槽侧视图
（ｂ）右挡板挡板传感器布置
Ｆｉｇ．２
图２水柱倒塌压力试验装置 Ｔｅｓｔｉｎｇ ａｐｐａｒａｔｕｓ ｆｏｒ ｗａｔｅｒ ｃｏｌｕｍｎ ｃｏＨａｐｓｅ
ｚ、Ｒ出
∞ 柏∞ 舳∞ 加０
０．０ ０．１ ０．２ ０．３ ０．４ ０．５ ０．６ ０．７ ０．８ ０．９ １．０
时间／ｓ （ａ）黏滞系数对压力的影响
舳∞ 加∞鲫∞ 加ｏ
Ｕ．０ Ｕ．１ ０．２ ０．ｊ ０．４ ０．５ ０．６ ０．７ Ｕ．８ Ｕ．９ １．Ｕ
时间／Ｓ
（ｂ）摩阻项对压力的影响
图５水柱倒塌后右挡板压力过程计算结果与
·１０３５·
由于各传感器采样时间间隔短暂，故挡板总压力误差 不大．而单宽压力则通过将总压力除以挡板宽度求得． ２组实验的挡板单宽压力时间过程分别如图４和图５ 所示．从２组实测结果来看，压力过程都呈现２个峰， 时间间隔约０．４５ Ｓ．水柱高度较小的第１组的第２个 峰较第１个峰明显要大，两峰相距时问间隔亦略短；第 ２组两峰压力值大致相同．
０．０ ０．１ ０．２ ０．３ ０．４ ０．５ ０．６ ０．７ ０．８ ０．９ １．０ 时间／Ｓ
图４水柱倒塌后右挡板压力过程计算结果与
实测结果对比（第１组次）
Ｆｉｇ．４
Ｃｏｍｐａｒｉｓｏｎ ｂｅｔｗｅｅｎ ｓｉｍｕｌａｔｅｄ ａｎｄ ｍｅａｓｕｒｅｄ ｒｅｓｕｌｔｓ ｏｆ ｐｒｅｓｓｕｒｅ ｐｒｏｃｅｓｓ ｏｎ ｒｉｇｈｔ ｂｏａｒｄ ａｆｔｅｒ ｗａｔｅｒ ｃｏｌｕｍｎ ｃｏｌｌａｐｓｅ（Ｃａｓｅ １）
Ｆ，＝一卷ｌ（口ｉ）Ｊ（鄱ｔ）
（２０）
式中：也是粒子半径；ｃＤ是摩阻系数，算例中取ｃＤ＝
１．５３［９１；（比ｉ）是该粒子处空间平均速度向量．
１．８数值计算过程 ＭＰｓ采用两步投影法，计算过程分显式和隐式两
阶段． 显式阶段：忽略动量方程（２）中的压力梯度项，仅
考虑体积力、黏滞力和摩阻力，对粒子的速度和坐标进 行第１次校正，速度增量为
数值实验结果表明，摩阻和黏滞项对压力计算结 果影响较大．以下结合实测结果分别予以讨论．
以第１组１８ ｃＩｎ水柱实验为例，将考虑和不考虑 摩阻２种情况的压力计算结果与实测结果进行了对比 （图４）．从计算结果可以看出，未加入摩阻项时，压力 峰尖锐，但峰点相位较实测有所提前．加入摩阻项后， 压力峰值略有减小，但峰点（Ａ点和Ｂ点）相位与实测 更接近．图６为峰点Ａ和曰时刻的瞬时流场分布图． 峰点Ａ对应流体与挡板前锋接触的时刻；峰点曰为流 体主体与挡板作用的时刻．
水柱倒塌后对挡板产生的总撞击力通过将各个传 感器压力乘以该传感器分担的挡板面积后相加获得，
万方数据
（ａ）水柱高度１８ ｃｍ
Ｆｉｇ．３
（ｂ）水柱高度２４ ｃｍ 图３粒子初始位置 Ｉｎｉｔｉａｌ ｃｏｎｆｉｇｕｒａｔｉｏｎ ｏｆ ｐａｒｔｉｃｌｅｓ
２００６年９月
李绍武等：ＭＰＳ水流数学模型在冲击压力计算中的应用
计算粒子数密度和压力梯度时采用较小的核子尺 寸ｒ。＝２．１ｄ。，可减少自由表面粒子数量；计算黏滞项 和压力二阶导数项时采用较大核子尺寸ｒ。＝３５ｄｏ，有 利于压力在较大范围进行分配．时间步长出与粒子尺 寸有关，粒子尺寸减小，时间步长应减小，根据Ｃｏｕｒａｎｔ 条件Ｈ ３计算中取Ａｔ＝０．０００ ４ Ｓ．粒子的初始压力由静 水压力给出． ３．２计算结果讨论
万方数据
万方数据
２００６年９月
李绍武等：ＭＰＳ水流数学模型在冲击压力计算中的应用
·１０３３·
计算各方向的粒子数密度和梯度，而不用空间维数来
补偿‘８｜，（即聊ｉ＝去荟［格（‰一
ｔ，ｔ）彬（Ｉ‘一‘ｆ）ｌ
（９）
式中
ｎｉ，。＝荟［！彳！；：｝三＿；孚伽（１ｌ一，；１）］
后＝菇，ｙ，ｚ
（１０）
１．４扩散模型
Ａ：坦塑
。㈣
上训（ｒ）出
由式（１１）和式（１２）得到
（嘞ｉ ２券荟∽圳训，：『一Ｉ） （１４）
用积分形式求出的Ａ和ｎ。Ａ是常数，未反映周围 粒子分布的不均匀差异．考虑到粒子分布的空间不均 匀性，Ａ离散计算形式为
∑伽（１
２
ｒＪ一‘１）Ｉ，：『一‘ｌ
卜—ｉ万Ｆ而一 ｑ５’
此外，‰亦用ｎｉ代替．结合ｎｉ的表达式（７），扩散 模型式（１４）随１可进一步写成
△砧ｉ＋＝△ｆ（Ｆ＋∥铲口？＋Ｆｒ） 粒子的速度和位置修正为
（２１）
谚＋１＝Ｈ，？＋÷（进＋１＋比？）△￡
（２３）
隐式阶段：经校正的粒子坐标使瞬时粒子数密度 ｎｉ＋偏离初始值凡。，通过求解压力泊松方程保证流体 的不可压缩性，同时求出新时刻的压力场【＿７｜，按式 （２４）求二次校正速度增量△即；一，按式（２５）和式（２６） 对粒子速度和坐标作二次校正，完成一个时间步．
取Ⅳ＝３．５０×１０“ｎ１２，考察了摩阻项对压力计算 结果的影响（图５（ｂ））．从计算结果看，尽管考虑摩阻 后峰值处略有损失，但总体结果与实测基本吻合．
万方数据