Page 78 - 水利学报2021年第52卷第6期

式中： q bk 为 k 粒径泥沙单宽推移质输沙率， m /s ；k 为粒径分组； u bk 为 k 粒径推移质运动速度，
2
m/s ； z 为河床高程，m； p ′ 为河床泥沙的孔隙率； α ′ ，α ′ 为推移质在水流和重力沿岸坡分量
b
m
by
bx
综合作用下的方向向量分量。根据不平衡输沙的经验公式，推移质的运动方程为：
æ q bk ö
∂ ç u ÷ ∂( α ′ q ) ∂( α ′ q )
è bk ø + bx bk + by bk + 1 （7）
∂t ∂x ∂y L bk (q - q *bk ) = 0
bk
式中： L 为 k粒径泥沙不平衡输移长度，m； q *bk 为 k粒径泥沙平衡输沙下的单宽推移质输沙率， m /s 。
2
bk
岸坡上的泥沙运动受重力和水流拖曳力的综合影响，推移质泥沙运动方向向量的修正计算公式为 ［15，72］ ：
α ′ bx = α + βsinφ x （8）
bx
α ′ by α by + βsinφ y
式中： α 、α by 为水流方向向量的分量； β 为考虑岸坡和床面泥沙性质的经验性修正系数，可以由
bx
实测资料率定或经验公式计算 ［73-76］ ； φ 、φ 分别为床面沿两个方向的倾角。
y
x
3.4 河冰动力学模块 河冰动力学模块的主要控制方程为河冰质量守恒方程、河冰面密度连续性方
程和动量方程 ［64］ 。河冰厚度和面密度由河冰连续性方程计算，河冰运动速度由动量方程计算。河冰
的动量方程、质量守恒方程、河冰质量及河冰面密度的连续性方程分别为：
      




M dV = R + F + F + G （9）
dt a w

dM + MÑV = E （10）
dt m
M = ρ Ct i （11）
i
dC  （12）
dt + CÑV + R - E = 0
a
a


   
式中：M 为单位面积冰的质量，kg； V 为冰速向量， m/s ； R 为冰内部的相互作用力，N； F 为风
a

  

作用在冰上的力，N； F 为水流作用在冰上的力，N； G 为重力在水平方向的分力，N； E 为单位
m
w
面积冰的增长速率，由热力生长消融或外部源汇控制，kg/s； ρ 为冰的密度， kg/m ； t 为冰的厚
3
i
i
-1 -1
度， m ； R 为机械作用引起的冰面积变化， s ； E 为冰面积的热力生长消融速率， s 。为了计
a
a
算河冰之间的内应力，河冰模块采用黏弹塑性应力应变方程计算冰块的内应力，采用摩尔-库伦屈服
准则计算冰塞冰坝中冰块的屈服及冰与河岸的临界摩擦力 ［77］ 。
式（9）—式（12）假设河冰颗粒为连续介质，采用基于拉格朗日场的光滑粒子法求解方程组。光滑
粒子法不受网格的限制，能准确计算冰盖前沿河冰下潜、冰塞体中冰块的坍塌变形及冰塞的冲蚀和
释放等复杂过程，避免了复杂水冰作用下离散格式的不稳定。
3.5 岸滩侵蚀模块 自然河流会出现周期性的岸坡冲刷变陡、河岸失稳破坏、坍塌土体落淤再平
衡的演变过程。大部分研究采用泥沙休止角判断非黏性沙岸坡的稳定与否，通过原型观测或参数
率定确定岸坡稳定的临界坡角，即泥沙休止角 ［78-80］ 。Zech 等针对淹没和非淹没岸坡不同含水量的岸
坡稳定性，提出双泥沙休止角法：采用稍大的动泥沙休止角判断岸坡的失稳与否，采用静泥沙休
止角决定崩岸坡面再平衡的稳定坡面，采用两组不同的动、静泥沙休止角分别计算水面上下的崩岸
过程 ［81］ 。本文岸滩侵蚀模块采用双泥沙休止角法计算冰盖及流凌刮擦影响下的岸滩侵蚀、崩塌和再
平衡过程 ［71］ 。计算方法如下。
（1）临界失稳点判断。如图 3 计算岸坡各个网格点的坡角，判断坡角是否大于相应的动泥沙休止
角，并确定临界崩塌破坏点的位置。采用式（13）计算坡面上各点的坡角：
-1æ z - z ö
α = tan ç bi + 1 bi ÷ （13）
i
è y i + 1 - y i ø
式中 α 为横断面 i 点的坡角。通过对比各点坡角与动泥沙休止角的大小，确定临界失稳点位置，该
i
点以上为不稳定坡面，以下为稳定坡面。图 3 中 N 为横断面网格点总数， N 为临界失稳点。
c
T
— 706 —