一种结冰水库回水区冰塞反演及预测方法与流程

1.本发明属于水利水电技术领域，涉及结冰水库回水区冰塞研究，尤其涉及一种基于水
‑
冰动力学模型耦合的结冰水库回水区冰塞反演及预测方法。


背景技术：

2.寒区河流上水库的修建，使库区及上下游河道的水力学特征发生显著变化，并直接影响到原河道的冰情过程。特别是水库回水末端，其水力学条件由河流过渡到水库，水面比降减小，水流流速降低，入库冰凌在此堆积；尤其是在结冰期容易产生冰塞，这将造成回水末端及部分上游河道水位壅高，严重时淹没两岸农田、房屋，对人民生命和财产安全产生直接威胁。水库库尾发生冰塞事件后，通常情况下采取低水位运行方式来减轻或消除灾害影响，其作用效果视具体情况而定，但较大程度直接影响到水库的发电效益。
3.尽管水库回水区冰塞现象普遍存在，其与河道冰塞在形成和演变机理方面具有一定的相似性，但影响要素更多、边界条件更为复杂而有所差别，具有河冰和水库冰的双重属性，采用河流冰塞理论难以直接描述水库回水区冰动力学变化和演变规律。限于水库回水区冰凌过程的复杂性和水库冰的研究进展，目前主要采用经验公式方法或静态冰坝理论对冰塞堆积体大小及水位壅高幅度进行估算，然而这些方法忽略了冰塞的动力过程，难以确定冰塞发生的时间和地点，同时更无法反映水库调度方式对回水区冰塞发展过程的影响。由于缺乏对水库回水区冰塞成因及过程的深入认识和系统有效的防凌调度方案，现阶段仍需要花费人力和物力进行现场实时追踪监测，以应对复杂的凌情变化并及时作出合理的防凌减灾措施。
4.因此，针对寒区工程的实际需求，开展水库回水区冰塞动力学的作用机理和演变规律研究，可为寒区水库工程设计、冰期调度运行和防凌减灾等方面提供重要的科技支撑，而建立能够同时考虑冰
‑
水动力耦合的模型，实现对结冰水库回水区冰塞反演和预测，是研究和解决这一问题急需有效的必要手段之一。


技术实现要素：

5.本发明的目的是为了克服现有技术中的不足，提供一种结冰水库回水区冰塞反演及预测方法，基于所建立考虑冰影响下的水动力和热力过程、表面冰动力、颗粒冰堆积层输移与堆积、堆冰头部位置的水深平均的平面二维水
‑
冰动力耦合模型，可以反演、预测和评估结冰水库回水区冰塞体及其引起的水位壅高的时空变化过程，以及水库回水区冰塞壅水对河道的影响程度和影响范围。
6.为达到上述目的，本发明采取以下技术来实现。
7.本发明提供的结冰水库回水区冰塞反演及预测方法，以结冰水库水域为研究对象，将其沿水流方向和宽度方向划分为若干三角形单元体，按照以下步骤对结冰水库回水区冰塞反演及预测：
8.s1给定水库边界条件以及水库初始水力、水温分布、冰场和堆冰头部位置条件；
9.s2利用水动力模块和热力模块获取当前时刻的水库回水区水力、水温以及水内悬浮颗粒冰浓度分布；
10.s3遍历水库，识别流速和水温与前一时刻相比，是否发生变化，若已发生变化，更新岸冰边界，然后进入步骤s4；否则直接进入步骤s4；
11.s4利用表面冰模块获取当前时刻的水库回水区表面冰压力、表面冰浓度和表面冰厚度分布；
12.s5遍历水库，判断水库表面各单元体节点的表面冰浓度是否大于设定的最大冰浓度，若是，表明表面冰发生机械增厚，修正表面冰厚度，然后进入步骤s6；若不是，则直接进入步骤s6；
13.s6遍历水库，利用颗粒冰堆积层输冰能力模型对存在表面冰的单元体节点下水流强度所能支撑的颗粒冰堆积层输移能力进行计算；
14.s7遍历水库，利用颗粒冰堆积层输移模型以及颗粒冰堆积层厚度变化模型，分别获得当前时刻的颗粒冰堆积层的冰流量和颗粒冰堆积层厚度分布；
15.s8在当前时刻基础上，增加时间步长；若增加后的时间不大于设定时间，返回步骤s2；若增加后的时间大于设定时间，当前程序结束，得到水库的水力、水温以及冰场分布情况。
16.本发明提供的结冰水库回水区冰塞反演及预测方法，是基于建立的水
‑
冰动力耦合模型实现的。所述水
‑
冰动力耦合模型包括水动力模块、热力模块、表面冰模块以及颗粒冰堆积输移模块。
17.(1)水动力模块
18.所述水动力模块包括建立沿水库水深平均的平面二维连续性模型和动量模型。
19.所述连续性模型为：
[0020][0021]
所述动量模型为：
[0022][0023][0024]
式中：x、y和t分别为空间和时间变量；q
x
、q
y
分别为单宽流量沿x、y方向的分量；ε
ij
是涡粘系数；τ
s
、τ
b
分别为水流与冰面和水流与河床的切应力；h为总水深；h
t
为表面冰或颗粒冰堆积层下的水深；η为水面高程。t
i
′
ce
为淹没冰厚，t
i
′
ce
＝ρ
ice
/ρ
·
t
ice
，t
ice
为表面冰厚度t
s
与颗粒冰堆积层厚度h
f
之和，ρ、ρ
ice
分别为水和冰的密度；n为表面冰浓度。
[0025]
(2)热力模块
[0026]
所述热力模块包括水温模型、水内悬浮颗粒冰模型。
[0027]
所述水温模型为：
[0028][0029]
所述水内悬浮颗粒冰模型为：
[0030][0031]
式中，t
w
为水温；c
v
为水内悬浮颗粒冰浓度；c
p
为水的比热容；ν
t
为水内悬浮颗粒冰离散系数；φ为水表面热通量；λ为水内悬浮颗粒冰与颗粒冰堆积层发生交换的体积率，
[0032]
为由于水体过冷却产生的水内悬浮颗粒冰浓度变化，可用以下公式确定：
[0033][0034]
式中：k
w
为水的热传导系数；为单个水内悬浮颗粒冰颗粒与水体的热交换nusselt数；n
f
为努赛尔数，n
f
＝c
v
/v0，为单位体积内水内悬浮颗粒冰的晶体个数；v0为水内悬浮颗粒冰颗粒的平均体积，α
o
＝πd
f
·
d
e
，为正对α轴的单个水内悬浮颗粒冰颗粒表面积，d
f
为水内悬浮颗粒冰颗粒α轴的长度，d
e
为水内悬浮颗粒冰颗粒的厚度；e
c
为冰晶单位面积与水的热交换率；l
ice
为水内悬浮颗粒冰的增长所释放的潜热；s
f
为由于水内悬浮颗粒冰增长引起的水温变化项，可简单表示为：
[0035]
s
f
＝
‑
p
t
t
w
ꢀꢀ
(1.5)
[0036]
式中，c
p
为水的比热容。
[0037]
在一定水流条件下，水内悬浮颗粒冰可能发生上浮并沉积在表面冰或颗粒冰堆积层底部，同时颗粒冰堆积层可能受水流剪切力的作用发生侵蚀，因此水内悬浮颗粒冰的浓度将发生变化，其质量守恒可用下式表示：
[0038][0039]
式中，式中，为到达表面冰或颗粒冰堆积层底部的水内悬浮颗粒冰的沉积概率；为单位面积上颗粒冰堆积层被侵蚀的定量系数；v
b
为水中悬浮颗粒冰的上浮速度；h为表面冰的等效厚度，h＝t
s
(1
‑
p)h
f
，p为颗粒冰堆积层孔隙率。
[0040]
在河流或水库中，还经常存在岸冰形态，它们的生成与当地的气象、水力条件密切相关，也是重要的冰边界。目前还没有成熟的数学模型来确定静态岸冰的生成，主要采用的是经验模式来确定。
[0041]
matousek根据ohre等河流现场观测总结得到河流表层水温的经验关系式：
[0042]
[0043]
式中：t
w,s
是引入的水表面温度；t
w
为水温，由公式(1.3a)计算得到；φ
wa
为水体向大气的失热量；为风速；b为经验系数，与风吹方向的水面宽度b相关；为水流速度。
[0044]
基于matousek提出的表层初始薄冰形成与飘移的水力学经验判别条件，提出水库静态表层初始薄冰盖形成条件为：
[0045][0046]
式中，v
cr
为水流临界速度，t
w,s
为水温临界值。
[0047]
(3)表面冰模块
[0048]
表面冰的运动主要受到水流、风力、重力、冰块间作用力以及它们与河岸间的相互作用力。本发明提供的表面冰动力学控制的表面冰模块包括运动模型、质量守恒模型、冰面积守恒模型和vep模型。
[0049]
所述运动模型、质量守恒模型和冰面积守恒模型为：
[0050][0051][0052][0053]
式中：为冰块加速度；m
ice
为单位面积冰质量，m
ice
＝ρ
ice
nt
ice
；ρ
ice
、n和t
s
分别为表面冰密度，表面冰浓度和表面冰厚；为冰块之间的作用力；为表面冰冰速；为风对表面冰的拖曳力；为水流拖曳力；为重力沿水流方向分力，g
x
、g
y
分别为重力在x、y方向的分量，r
a
为由于机械再分配而引起的冰面积浓度变化。
[0054]
可表达为：
[0055][0056]
式中：σ
xx
、σ
yy
分别为冰块之间的正应力；σ
xy
、σ
yx
分别为冰块之间的剪应力。
[0057]
冰的内部作用力需要构建一个本构模型来描述。因此，一种广泛的粘弹塑性模型(vep)被引入用来计算冰的内部应力。在二维模式下，vep模型可以写为：
[0058][0059]
式中：k、g分别为河冰的体积弹性模量、剪切弹性模量；ξ
v
、η
v
分别为河冰的体积粘性系数、剪切粘性系数；p为压力项；δ
ij
为克罗内克尔符号；k＝x，y；
u、v分别为表面冰速度沿x、y方向的分量。
[0060][0061][0062]
式中，γ为泊松比，取0.3；；e为杨氏模量，取值为1.0
×
109pa。
[0063]
p为压力项，通过扩展静态冰坝的构成基本关系，采用以下表达式来描述：
[0064][0065]
式中：θ为冰的内摩擦角，θ＝46
°
；n
max
为冰的最大聚集程度，n
max
＝0.6；λ为经验参数，λ＝15。
[0066]
冰的当前含量极限不能超过最大值n
max
，当达到最大值时，就会发生机械增厚，冰的厚度将修正为：
[0067][0068]
式中的n＞n
max
。
[0069]
(4)颗粒冰堆积输移模块
[0070]
所述颗粒冰堆积输移模块包括输冰能力模型、颗粒冰堆积层输移模型和颗粒冰堆积层厚度变化模型。
[0071]
水流中悬浮颗粒冰，包括由于水流中紊动、冰块之间碰撞造成的表面冰脱落形成和水流自身由于失热、过冷却形成。这些颗粒冰在表面冰或已形成的颗粒冰堆积体底部发生沉积，进而发生输移。表面冰堆积体下颗粒冰堆积层的输移能力，与当地水流条件和冰粒特征有关，水流强度所能支撑的颗粒冰堆积层输移能力(也为水流的输冰能力)采用以下公式确定：
[0072][0073]
式中：q
uc
为单位宽度内颗粒冰堆积层的冰流量通过能力；θ无因次水流强度；θ
c
为无因次临界水流强度，θ
c
＝0.041；d
n
为冰颗粒标称尺寸，取值为0.0001m；s为冰颗粒与水的密度比，s＝ρ
ice
/ρ；f为颗粒冰上浮速度系数，对于球状颗粒近似等于1.0。
[0074]
θ可由以下表达式确定：
[0075][0076]
式中：u
*ice
为表面冰或颗粒冰堆积层下侧的剪切速度。
[0077]
颗粒冰堆积层输移模型的二维形式为：
[0078]
[0079]
式中：q
u
为单位宽度内颗粒冰堆积层的冰流量；u
f
为冰盖下颗粒冰堆积层流速，根据van rijn(1984b)和wu等(2006)所定义的公式确定；l
x
为颗粒冰堆积层的自适应长度；q
ux
为x方向单位宽度颗粒冰堆积层的冰流量，q
ux
＝q
u
cosθ；q
uy
为y方向单位宽度颗粒冰堆积层的冰流量，q
uy
＝q
u
sinθ；θ为流动方向与x正方向的角度；
[0080]
颗粒冰堆积层厚度h
f
变化模型的方程可表达为：
[0081][0082]
q
f
为水内悬浮颗粒冰层与表面冰底部或颗粒冰堆积层的单位面积交换率，可表示为：
[0083][0084]
式中：p为颗粒冰堆积层堆积的孔隙率；h
f
为颗粒冰堆积层堆积厚度；n
ice
为颗粒冰堆积层的冰面积浓度，n
ice
＝q
u
/(h
f
u
f
)。
[0085]
上述结冰水库回水区冰塞反演及预测方法，以结冰水库水域为研究对象，将其沿规定的直角坐标系划分为若干三角形单元体，以每个单元体所在节点的水力、水温、冰场表征其所在位置的水力、水温和冰场情况。所述水力分布情况包括结冰水库的流量、水深、流速随时间变化过程；所述水温分布情况包括结冰水库的库区水温随时间变化过程；所述冰场分布情况包括岸冰、水内悬浮颗粒冰、颗粒冰堆积层的时空分布。对水动力学模块、热力模块和颗粒冰堆积输移模块进行处理的方式为，采用有限元法对模型进行离散求解处理，对表面冰模块进行的处理方式为，采用sph(smoothed particle hydrodynamics，光滑粒子流体动力学)方法进行求解处理。在给定时间段内，按照上述步骤s1
‑
s12，便可得到计算区域内水力(水深、流量、流速)、水温、冰情(岸冰、水内悬浮颗粒冰浓度、表面冰浓度及厚度、颗粒冰堆积层冰流量及厚度)分布情况。
[0086]
上述步骤s1中，所述水库边界条件可以通过本领域常规手段来确定，其包括上游入流边界条件、下游出流边界条件、表面边界条件、河床底部边界条件、河岸边界条件：
[0087]
101、上游入流边界条件
[0088]
上游水流为全断面入流，可采用流量边界条件，也可采用水位边界条件；上游边界的水温采用实测或给定的水温，水内悬浮颗粒冰浓度采用实测或给定的浓度值，表面冰厚度、浓度分别采用实测或给定的厚度、浓度值，表面冰的速度设定与水流速度一致。
[0089]
由于采用sph方法计算表面冰，需要在上游入流边界断面处释放冰粒子，可在选取的断面输冰宽度内采用总量控制的随机方法给定。
[0090]
102、下游出流边界条件
[0091]
下游水流边界为全断面出流，下游边界条件可选择给定水位边界条件(对应上游边界为流量边界条件)，也可选择给定流量边界条件(对应上游边界为水位边界条件)。
[0092]
水温、水内悬浮颗粒冰、颗粒冰堆积层出口边界的梯度为零。
[0093]
103、水库表面边界条件
[0094]
水流表面的水动边界条件，满足自由表面运动学和动力学边界条件。当水流表面存在流冰时，还应考虑冰对水流的阻力影响。
[0095]
此外，水面不存在冰时，主要考虑气
‑
水的热交换过程。水面若有冰，主要考虑冰
‑
水的热交换过程。
[0096]
104、水库河床底部边界条件
[0097]
库底采用无滑移边界条件，且为绝热边界。
[0098]
105、河岸/岸冰边界
[0099]
当冰向边界移动时，河床/岸冰边界对流冰的摩擦阻力需要考虑，其计算遵循动态库仑准则：
[0100]
f
f
＝f
c
f
n tanφ
b
[0101]
式中，f
f
为冰与边界之间的摩擦力；f
c
为黏合力，假设为0；tanφ
b
为动力学摩擦系数；f
n
为作用于冰边界上的垂直分量，该数值可通过计算冰内部垂直于边界的力来获得。
[0102]
上游流冰堆积至回水区，表面冰堆冰头部所在位置断面的平均水深与流量满足公式(1.24)的关系。
[0103]
h
′
＝0.0134q 1.9
ꢀꢀ
(1.24)式中，h
′
为表面冰堆冰头部位置所在断面的平均水深，q为断面流量。
[0104]
本发明中，以表面冰堆冰头部所在位置作为结冰水库回水区冰塞反演或预测过程中的约束条件。
[0105]
所述水库水力、水温分布初始条件，通过上游边界采用不变的初始流量和水温、下游边界采用水动力模块(公式1.1
‑
1.2b)和热力模块(公式1.3)从参考水位(网格划分时给定)变化至给定初始水位的计算来获取拟定初始条件下单元体节点的水深、流速、水温等相关信息。
[0106]
上述步骤s2中，依据前一时刻的水库水力分布和冰场分布作为当前时刻水库水力分布和冰场分布的初始值，依据水动力模块的连续性模型、动量模型以及热力模块的水温模型、水内悬浮颗粒冰模型，然后按照以下分步骤得到当前时刻的水库水力参数、水温分布和悬浮颗粒冰浓度：
[0107]
s201通过连续性模型(公式1.1)和动量模型(公式1.2a
‑
1.2b)获取水库回水区的水深、流量、流速参数；
[0108]
s202通过水温模型(公式1.3a)获取水库回水区的水温分布；
[0109]
s203通过水内悬浮颗粒冰模型(公式1.3b
‑
1.6联立)获取水库回水区的水内悬浮颗粒冰浓度分布。
[0110]
上述步骤s3中，遍历水库表层每个单元体，将每个单元体各节点的流速和水温与前一时刻相比，从而识别每个单元体各节点的流速和水温是否发生变化。当t
w
>0℃，t
w,s
≤t
cr
，且时，表明岸冰生成，该位置不能让表面冰通过。此时将满足条件的节点构成的单元体设定为岸冰，从而更新岸冰边界。
[0111]
上述步骤s4中，采用sph方法求解表面冰模块包含的运动模型、质量守恒模型、冰面积守恒模型和vep模型(公式1.9
‑
1.16)得到当前时刻的水库回水区表面冰压力、表面冰浓度和表面冰厚度分布。
[0112]
上述步骤s5中，各单元体节点当前冰含量极限不能超过冰的最大聚集程度(即最大值n
max
)，当达到最大值时，就会发生机械增厚，冰的厚度将修正为：
[0113]
[0114]
式中，这里的n＞n
max
。
[0115]
再基于修正后的冰厚度，执行步骤s6。若单元体节点当前冰含量没有超过冰的最大聚集程度，可以直接执行步骤s6。
[0116]
上述步骤s6中，对于存在表面冰的单元体节点(即单元体节点表面冰浓度大于0)，按照公式(1.18)计算相应单元体节点下方颗粒冰堆积层的输冰能力。
[0117]
上述步骤s7中，遍历水库，利用颗粒冰堆积层输移模型(公式1.20)和颗粒冰堆积层厚度变化模型(1.21)得到各单元体节点颗粒冰堆积层的冰流量和颗粒冰堆积层厚度分布。
[0118]
上述步骤s8中，在当前时刻基础上增加一个时间步长δt，并判断增加后的时间是否达到设定的时间段上限t
end
。若没有达到设定的时间段上限t
end
，返回步骤s2，将利用当前时刻得到的水力、水温分布、冰场各参数对各变量进行赋值，作为下一时刻的水力、水温、冰场参数变量初始值；重复上述步骤s2
‑
s8的操作；若达到设定的时间段上限，则程序结束。对计算得到的水力、水温、水内悬浮颗粒冰浓度、表面冰压力、表面冰浓度、表面冰厚度、颗粒冰堆积层冰流量、颗粒冰堆积层厚度数据进行汇总统计，便可得到水库的水力、水温和冰场分布情况。所述水力分布情况包括结冰水库的流量、水深、流速随时间变化过程；所述水温分布情况包括结冰水库的库区水温随时间变化过程；所述冰场分布情况包括岸冰、表面冰、颗粒冰堆积层、水内悬浮颗粒冰的时空分布。
[0119]
为了实现对结冰水库回水区冰塞的准确反演或预测，本发明对水库水力参数、水温、水内悬浮颗粒冰浓度和颗粒冰堆积层的求解采用的是一种改进的显式时间积分方法
‑
四阶龙格
‑
库塔法，其时间步长应满足：
[0120][0121]
其中，δl为单元体的最小长度；h
max
为最大水深。
[0122]
表面冰动力学的求解仍采用四阶龙格
‑
库塔法，冰块计算的时间步长由以下公式确定：
[0123][0124]
其中，a
k
、c
k
为由冰属性决定的加速度和声速；β是courant数，值为1.0。
[0125]
与现有技术相比，本发明具有以下有益效果：
[0126]
(1)本发明方法，在考虑水动力
‑
冰动力
‑
地形耦合作用下，对建立的水深平均的平面二维水
‑
冰动力耦合模型进行求解，得到结冰水库的水力分布、水温分布和冰场分布，从而实现对结冰水库库尾回水区的冰塞进行反演或预测；由于考虑了水动力和热力模型、表面冰动力、水内悬浮颗粒冰、颗粒冰堆积层输移、堆冰头部位置等影响因素，可以更加准确地反演或预测不同水库运行方式下、不同来冰条件下结冰水库库尾回水区的冰塞体形态、大小与引起的壅水水位，并可获取动力演进的时空变化过程。
[0127]
(2)本发明考虑了水库回水区冰塞体在水动力、表面冰动力、颗粒冰堆积层堆积与输移、静态岸冰等不同类型冰水之间的动力学及其与地形等条件耦合情形，能反映回水区冰塞体冰水之间的质量、动量守恒及传递过程。
[0128]
(3)本发明将水库回水区冰塞体的来源组成分为表面冰和颗粒冰堆积层，颗粒冰
堆积层对冰塞体冰量的贡献与表面冰相当；颗粒冰堆积层的冰源，不仅仅考虑了由于水流流动过程中由于失热和进一步过冷却产生的水内冰，还包括由于表面冰运动过程中受水流紊动和冰团之间碰撞脱落后回到水中二次处于悬浮的颗粒冰。
[0129]
(4)表面冰堆冰头部位置对于水库库尾回水区冰塞的演变起着至关重要的作用，目前并没有完善的基础理论和数学模型对其进行确定，传统方法取断面平均流速0.3～0.4m/s所处的位置；本发明方法根据万家寨水库现场原位监测资料，建立了水库回水区表面冰堆冰头部位置与断面水深、流量的经验关系，该经验关系式更能准确表达水库回水区表面冰堆冰头部在地形和水库调度方式下的位置。
附图说明
[0130]
图1为结冰水库回水区表面冰、颗粒冰堆积层与水内悬浮颗粒冰的分布示意图。
[0131]
图2为本发明采用的水
‑
冰动力耦合模型框架图。
[0132]
图3为本发明提供的结冰水库回水区冰塞预测及反演方法流程示意图。
[0133]
图4为实施例中万家寨水库回水区水位计算值与观测值的对比图。
[0134]
图5为实施例中万家寨水库回水区冰塞体与水位随时间的演变示意图；其中，(a)对应第7小时，(b)对应第27小时，(c)对应第46小时，(d)对应第63小时，(e)对应第80小时。
[0135]
图6为实施例中万家寨水库回水区形成的冰塞体稳定后的各断面冰厚计算值与观测值的对比图。
具体实施方式
[0136]
以下将结合附图对本发明实施例的技术方案进行清晰、完整的描述，显然，所描述实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施例，都属于本发明。
[0137]
本发明中所述的结冰水库回水区的研究范围指的是水库大坝至回水到天然状态的河段。
[0138]
实施例
[0139]
本实施例运用本发明所提供的结冰水库回水区冰塞反演及预测方法获取万家寨水库在2014
‑
2015年度冬季流凌前和流凌封河期的水温分布情况以及冰场分布情况。如图1所示，水库最表面为表面冰层分布区域，表面冰层下方为颗粒冰堆积层，水流中分布着大量的悬浮颗粒冰。进入表面冰堆积河段后，水内悬浮颗粒冰在短距离内上升到达表面冰底面形成颗粒冰堆积层，颗粒冰堆积层整体以与推移质泥沙相似的盖移质形式输移。颗粒冰堆积层的冰流量受水流的输冰能力影响，使表面冰底部发生堆积或侵蚀。当颗粒冰堆积层的冰流量大于水流的输冰能力时，表面冰底部就会发生堆积；当颗粒冰堆积层的冰流量小于水流的输冰能力时，表面冰底部就会发生冲刷侵蚀；当颗粒冰堆积层的冰流量等于水流条件所能支撑的输冰能力时，表面冰底部颗粒冰的沉积与冲刷达到平衡。
[0140]
首先采用三角形网格，在水库三维地形基础上进行网格划分。三角形网格的边大小可以根据水库水面宽进行设定，一般要求在最小河段宽度上需要分布8个以上的三角形网格。
[0141]
本实施例中在建立的x
‑
o
‑
y空间坐标中，将万家寨水库离散成由19395个节点(m，n)组成的35250个三角形网格，网格尺寸平均为35m，以每个节点的水温、冰场情况表征其所在位置的水温和冰场情况，m表示水流方向的位置坐标，n表示沿宽度方向的位置坐标。
[0142]
本实施例基于的水
‑
冰动力耦合模型如图2所示，其包括水动力模块、热力模块、表面冰模块以及颗粒冰堆积输移模块。水动力模块包括连续性模型(公式(1.1))和动量模型(公式(1.2a)
‑
(1.2b))。热力模块包括热力模型(公式(1.3a))、水内悬浮颗粒冰模型(公式(1.3b))。表面冰模块包括运动模型(公式(1.9))、质量守恒模型(公式(1.10))、冰面积守恒模型(公式(1.11))和vep模型(公式(1.13))。颗粒冰堆积输移模块包括颗粒冰堆积层输移能力模型(公式(1.18))、颗粒冰堆积层输移模型(公式(1.20))和颗粒冰堆积层厚度变化模型(公式(1.21))。
[0143]
上述各模型中的相关参数变量取值见表1所示。
[0144]
表1各模型中的相关参数变量
[0145][0146]
本实施例采用有限元法对水动力学模块、热力模块和颗粒冰堆积输移模块进行离散，得到各模型的离散方程。对表面冰模块的进行的处理方式为，采用sph(smoothed particle hydrodynamics，光滑粒子流体动力学)方法进行求解处理。
[0147]
如图3所示，本实施例按照以下步骤对万家寨水库的水温分布情况以及冰场分布情况进行反演：
[0148]
s1给定水库边界条件以及水库初始水力、水温分布、冰场和堆冰头部位置条件。
[0149]
本实施例万家寨水库所涉及的边界条件以及冰场的确定方法如前给出的101
‑
105所述。
[0150]
水库水力、水温分布初始条件，通过上游边界采用不变的初始流量和水温、下游边界采用水动力模块(公式1.1
‑
1.2b)从参考水位984.0m变化至给定初始水位(2014.11.05实测水位：976.1m)的计算来获取拟定初始条件下单元体节点的水深、流量等相关信息，水温分布采用水库沿程不同断面的实测点内插得到。
[0151]
水库水内悬浮颗粒冰浓度、颗粒冰堆积层的浓度及厚度、表面冰浓度及厚度的初始值均设定为0。
[0152]
表面冰堆冰头部所在位置根据公式(1.24)计算得到。
[0153]
s2利用水动力模块和热力模块获取当前时刻的水库回水区水力、水温以及水内悬浮颗粒冰浓度分布。
[0154]
本步骤依据前一时刻的水库水力参数、水温分布、冰场分布(包括水内悬浮颗粒冰浓度、表面冰压力、表面冰浓度、表面冰厚度、颗粒冰堆积层冰流量、颗粒冰堆积层厚度)各变量作为当前时刻的初始值(对于初始时刻，当前时刻水库水力参数、水温分布、冰场分布对应的是初始条件)，然后按照以下分步骤得到当前时刻的结冰水库水力参数、水温分布和水内悬浮颗粒冰浓度：
[0155]
s201通过连续性模型和动量模型获取水库的水力参数。
[0156]
本步骤中，依据边界条件和当前时刻结冰水库水力各变量初始值、水温分布相应变量初始值，通过有限元法对连续性模型(公式(1.1))和动量模型(公式(1.2a)和(1.2b))进行离散求解，得到当前时刻水流的总水深h(m,n)、表面冰或颗粒冰堆积层下的水深h
t
(m,n)和流量分布q
x
(m,n)、q
y
(m,n)。然后根据或可得到水流速度
[0157]
s202通过水温模型获取水库的水温分布。
[0158]
本步骤，依据边界条件、当前时刻的水温初始值以及步骤s201得到的水力参数对水温模型(式(1.3a))的离散方程进行求解，得到水温场参数，t
w
(m,n)。
[0159]
s203通过水内悬浮颗粒冰模型获取水库的水内悬浮颗粒冰浓度分布。
[0160]
本步骤中，依据边界条件、当前时刻水内悬浮颗粒冰浓度初始值以及步骤s201得到的水力参数和步骤s202得到的水温参数对水内悬浮颗粒冰模型(式(1.3b)
‑
(1.6)联立)的离散方程进行求解，得到当前时刻的水内悬浮颗粒冰浓度参数，c
v
(m,n)。
[0161]
s3遍历水库，识别流速和水温与前一时刻相比，是否发生变化，若已发生变化，更新岸冰边界，然后进入步骤s4；否则直接进入步骤s4。
[0162]
本步骤中，遍历水库表层每个单元体，将每个单元体各节点的流速和水温与前一时刻相比，从而识别每个单元体各节点的流速和水温是否发生变化。当t
w
>0℃，t
w,s
≤t
cr
，且且时，表明岸冰生成，该位置不能让表面冰通过。此时将满足条件的节点构成的单元体设定为岸冰，从而更新岸冰边界。
[0163]
s4利用表面冰模块获取当前时刻的水库回水区表面冰压力、表面冰浓度和表面冰厚度分布。
[0164]
本步骤中，采用sph方法，依据步骤s3更新后的岸冰边界，通过上述公式(1.9)
‑
(1.16)求解，便可得到当前时刻的水库表面冰压力p(m,n)、表面冰浓度n(m,n)和表面冰厚度t
s
(m,n)分布。
[0165]
s5遍历水库，判断水库表面各单元体节点的表面冰浓度是否大于设定的最大冰浓度，若是，表明表面冰发生机械增厚，修正表面冰厚度，然后进入步骤s6；若不是，则直接进入步骤s6。
[0166]
本步骤中，首先判断各单元体节点当前冰含量是否大于设定的最大冰浓度，若是按照公式(1.17)对表面冰厚度进行修改，再基于修正后的冰厚度，执行步骤s6；若不是，可以直接执行步骤s6。值时，就会发生机械增厚，冰的厚度将修正为：
[0167][0168]
式中，这里的n＞n
max
。
[0169]
若单元体节点当前冰含量没有超过冰的最大聚集程度，
[0170]
若表面冰浓度大于设定的最大冰浓度，表明表面冰发生机械增厚，(采用公式1.17)。
[0171]
若表面冰浓度不大于设定的最大冰浓度，表明表面冰没有发生机械增厚，直接进入步骤s6。
[0172]
s6遍历水库，利用颗粒冰堆积层输冰能力模型对存在表面冰的单元体节点下水流强度所能支撑的颗粒冰堆积层输移能力进行计算。
[0173]
本步骤中，对于存在表面冰的单元体节点，按照公式(1.18)计算相应单元体节点下方颗粒冰堆积层的冰输冰能力。若当地无因次水流强度超过了临界值θ
c
(＝0.041)，则计算得到当地水流强度所能支撑的颗粒冰堆积层的冰输冰量。若当地无因次水流强度小于临界值θ
c
(＝0.041)，则认为当地水流条件下可支撑的颗粒冰堆积层的冰输冰量为0。
[0174]
s7遍历水库，利用颗粒冰堆积层输移模型以及颗粒冰堆积层厚度变化模型，分别获得当前时刻的颗粒冰堆积层的冰流量和颗粒冰堆积层厚度分布。
[0175]
本步骤中，通过对颗粒冰堆积层输移模型(公式1.20)和颗粒冰堆积层厚度变化模型(1.21)的离散方程进行求解，得到当前时刻的得到各单元体节点颗粒冰堆积层的冰流量q
u
(m,n)和颗粒冰堆积层厚度h
f
(m,n)分布。值得注意的是，本实施例中设定表面冰底部不具备可侵蚀性。
[0176]
s8在当前时刻基础上，增加时间步长；若增加后的时间不大于设定时间，返回步骤s2；若增加后的时间大于设定时间，当前程序结束，得到水库的水力、水温以及冰场分布情况。
[0177]
在当前时刻基础上增加一个时间步长δt，并判断增加后的时间是否达到设定的时间段上限t
end
。若没有达到设定的时间段上限t
end
，返回步骤s2，将利用当前时刻得到的结冰水库水力参数、水温、冰场各参数对各变量进行赋值，作为下一时刻的水库水力参数、水温、冰场参数变量初始值；重复上述步骤s2
‑
s8的操作；若达到设定的时间段上限，则程序结束。对计算得到的水力、水温、水内悬浮颗粒冰浓度、表面冰压力、表面冰浓度、表面冰厚度、颗粒冰堆积层冰流量、颗粒冰堆积层厚度数据进行汇总统计，便可得到反演的结冰水库的水温分布情况以及冰场分布情况。
[0178]
图4给出了运用本发明提供的结冰水库回水区冰塞反演及预测方法对黄河万家寨水库反演得到的2014
‑
2015年度流凌封河前(畅流期)的模拟水位与实测水位的对比。坝前水位的计算值与观测值平均差异为
‑
0.01m，差异范围为
‑
0.28m～0.18m。水泥厂位于变动回水区，其计算值与观测值平均差异为
‑
0.06m。可以看出，本发明可以较好地模拟水库水位、入流流量同步变化条件下沿程水位随时间的变化过程，且具有较高的精度。
[0179]
图5给出了运用本发明提供的结冰水库回水区冰塞反演及预测方法对万家寨水库2014
‑
2015年度流凌封河期反演得到的库尾回水区冰塞体与水位随时间的演变过程图。第7小时(初始表面流冰已开始2小时)，库尾流冰基本到达回水区开始堆积，且该时段的初始冰盖(表面冰)已经形成。大约6小时后(距离初始流冰约8小时)，流冰开始到达水库回水区冰盖前沿并形成堆积，而库尾冰盖上游河道中的颗粒冰堆积层也同步到达表面冰下，由于流速和水流紊动的降低，颗粒冰堆积层开始上浮并累积在表面冰前沿的下游位置，表面冰下较小的流速致使冰盖下堆积的颗粒冰堆积层输运能力减缓，导致水内悬浮颗粒冰大量堆积。第27小时，回水区表面冰不断增厚，主要因为持续来冰的挤压作用，同时颗粒冰堆积层的堆积也随着冰塞的形成而向上游发展，冰塞头部因颗粒冰堆积层的大量堆积造成了当地水深的减小和流速的增加，由此增加了水对冰的拖曳力以及颗粒冰堆积层的输运能力。颗粒冰堆积层在表面冰下的堆积，以及表面冰前沿因挤压而发生的机械增厚，导致了上游水位壅高和冰输运能力的降低，加上弯道及河岸阻力作用，冰塞较快地向上游发展。第46小时，岔河口上游的冰塞延伸至水泥厂的情况，而岔河口下游的冰塞冰厚缓慢增长，这是因为其上游的冰塞大大减少了对下游冰的供应。由于水的拖曳力和冰的重力作用，此冰塞冰厚因冰内部的挤压而继续增厚。进一步，第63小时水泥厂水位到达了实测的最高水位。随后由于上游边界水位的持续降低(来水流量变小)，水泥厂和岔河口的水位也随着降低。第80小时，库尾河段的冰塞体和水位达到稳定状态，上游来冰不能继续进入该河段，以平铺方式向上游河段继续发展。
[0180]
图6给出了运用本发明提供的结冰水库回水区冰塞反演及预测方法对万家寨水库2014
‑
2015年度流凌封河期反演得到的库尾回水区不同位置冰塞体厚度(表面冰厚度与颗粒冰堆积层厚度之和)与稳封期观测值的对比。反演结果表明，断面wd54
‑
56之间，最大冰厚反演计算值比观测值略小；断面wd56
‑
60之间，平均冰厚模拟计算值比观测值略大；随着时间的推移，断面wd56
‑
60之间的颗粒冰堆积层逐渐向下游输移，从而使wd54
‑
56河段冰厚增加；断面wd61
‑
65之间，由于反演时段内的边界条件难以精确获取，导致模拟过程中冰塞体的发展与实际情况存在一定的差异。总体来看，库尾河段各断面最大冰厚、平均冰厚反演计算值与稳封期观测值在趋势上是一致的，且模拟计算最大冰厚位置与实测最大冰厚位置相同(在冰塞头部位置wd50)。
[0181]
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。