Page 17 - 水利学报2025年第56卷第4期

顿第二定律，颗粒运动的控制方程如下
ｄｕ
ｐ
＝ Ｆ ＋ Ｆ ＋ Ｆ ＋ Ｆ （３）
ｄｔ Ｄ Ｂ Ｐ Ｖ
式中：ｕ为颗粒速度；Ｆ 为曳力；Ｆ 为浮力；Ｆ 为压力梯度力；Ｆ 为虚拟质量力。
ｐ
Ｖ
Ｄ
Ｂ
Ｐ
颗粒受力的计算公式如下
３ ρ ｆ
Ｆ ＝－ Ｃ ｕ－ ｕ （ｕ－ ｕ） （４）
Ｄ Ｄ ｐ ｆ ｐ ｆ
ｐ
４ｄ ρ ｐ
－ ）
ｇ（ ρ ｐ ρ ｆ
Ｆ ＝ （５）
Ｂ
ρ ｐ
ρ ｆ
Ｆ ＝ ｕ!ｕ ｆ （６）
ｐ
Ｐ
ρ ｐ
ρ ｆ ｄ ｕ ｐ
Ｆ ＝ Ｃ Ｍ ( ｕ!ｕ－ ) （７）
Ｖ
ｐ
ｆ
ρ ｐ ｄｔ
为颗粒的密度；Ｃ 为曳力系数 ［１９］ ；Ｃ 为虚拟质量系数，取 ０．５；ｄ为颗粒粒径；ｇ为重力加
式中：ρ ｐ Ｄ Ｍ ｐ
速度。
本文采用 Ｇｒａｎｔ等 ［２０］ 提出的颗粒－ 壁面反弹模型模拟颗粒碰撞壁面前后的运动，如下
ｖ ２ ３
ｐ２
ｎ
ｅ＝ ＝ ０．９９３ － １．７６ θ ｒ ＋ １．５６ θ ｒ － ０．４９ θ ｒ
ｖ
ｐ１
（８）
ｕ ｐ２
３
２
ｅ ＝ ＝ ０ ．９８８ － １．６６ θ ｒ ＋ ２ ．１１ θ ｒ － ０ ．６７ θ ｒ
ｔ
ｕ
ｐ１
式中：ｅ为弹性恢复系数的法向分量；ｅ为切向分量；ｖ和 ｖ为碰撞前后颗粒速度的法向分量；ｕ 和
ｎ ｔ ｐ２ ｐ１ ｐ２
为颗粒对壁面的撞击夹角。
ｐ１
ｕ 分别为碰撞前后颗粒速度的切向分量；θ ｒ
［２１］
颗粒对壁面的磨损采用 ＧｅｎｅｒｉｃＥｒｏｓｉｏｎＭｏｄｅｌ 进行模拟，即
Ｃ（ｄ）ｆ（ α ）ｖ
Ｎ ｍ ｂ（ｖ）
ＥＲ ＝ ∑ ｎ Ｐ （９）
ｎ ＝１ Ａ ｆａｃｅ
式中：ＥＲ为磨损率，定义为颗粒对单位面积靶材造成的质量损失；Ｎ为颗粒撞击壁面的数目； ｍ 为
ｎ
－ ９
颗粒的质量流率；Ａ 为颗粒撞击网格单元的面积；Ｃ（ｄ）为颗粒粒径的函数，值为 １．８ × １０ ；ｂ（ｖ）为
Ｐ
ｆａｃｅ
颗粒相对速度 ｖ的函数，值为 ２．６；ｆ（ α ）为颗粒冲角的函数：
０．０４ α α ＜２０°
０．０８ ＋ ０．０２（ α － ２０） ２０° ≤α ＜３０°
α － ３０
ｆ（ α ） ＝ １ － ３０° ≤α ＜４５° （１０）
３０
０ ．１（ α － ４５ ）
０．５ － ４５° ≤α ＜９０°
４５
２．３ 边界条件及求解设置 本文采用标准离散格式离散压力，二阶迎风格式离散动量、湍流耗散率和比
耗散率等，求解方法采用 ＳＩＭＰＬＥＣ算法。由于 Ｊｈｉｍｒｕｋ水电站长期处于大开度的运行工况，因此本文选
取的计算工况为大开度（１８°）工况。进、出口均采用压力边界条件，其中进口为 ４０ｍ水头，出口为一个
３
标准大气压。根据水电站实测数据，计算采用的泥沙粒径为 ０．１ｍｍ，浓度为 ５００ｍｇ?Ｌ，密度为 ２６５０ｋｇ?ｍ。
２．４ 网格划分及网格无关性分析 如图 ２所示，本文采用 ＡＮＳＹＳＩＣＥＭ对计算域进行六面体结构化网
格划分。其中，由于转轮叶片区域几何外形的复杂性，采用 ＴＵＲＢＯＧＲＩＤ进行六面体非结构化网格划
分。对近壁处的网格进行加密处理，叶片部分的 ｙ ＋ 在 ３０以内，其余部分 ｙ ＋ 在 １００以内。
网格数量会影响计算精度和效率，需要进行网格无关性分析 ［２２ － ２３］ 。由图 ３（ａ）可知，随着网格数
目的增加，水轮机的效率和平均磨损率的变化渐趋稳定，考虑到计算精确性和成本，最终选择网格总
— ４ ３ ５ —