水库异重流发生及形成的特定位置，一般称为潜入点。异重流潜入点是判断异重流与明流流态的重要节点，其流速分布和含沙量分布的具体形式，与其上、下游断面的流速和含沙量的分布形式不同，属于过渡状态，此乃判别潜入点的重要依据，也是从前人研究获得的认识。流速分布是水流阻力状况的反映，或者说是水流能量消耗的反映，研究水库异重流潜入点垂线流速分布理论公式，进一步改进潜入条件，对水库调水调沙预案编制、水库优化调度及多沙河流水库规划设计等方面具有重要意义。

在研究水库异重流的流速与含沙量沿垂向分布规律方面，文献[1-2]给出了层流异重流的流速分布为抛物线型和紊流异重流流速分布形式，但异重流潜入点处流速复杂多变，其垂线分布形式尚未给出。围绕米德尔湖、官厅水库、刘家峡水库等水库收集的实测资料，中国科学院水利水电科学研究院^[3]结合水槽试验研究了不同进口含沙量的异重流潜入条件表达式。Savage和Brimberg^[4]通过理论分析底坡和交界面与河底阻力的影响，得到修正弗劳德数。Farrell和Stefan^[5]采用水槽试验和数学模型研究异重流的潜入和运动。张俊华^[6]、李书霞、夏军强^[7]、李涛^[8]等采用实测资料分析、水槽试验、物理模型实验，分别研究了异重流潜入点条件、潜入条件中的重力修正系数随含沙量变化等问题。也有学者从异重流潜入物理机制出发研究异重流潜入理论模型，赵琴等^[9]利用动量方程和能量方程得到了浑水异重流潜入理论模型。李书霞、夏军强等^[7]分析了异重流的动量修正系数与异重流潜入处的流速分布关系，在水库异重流潜入条件中假定动量修正系数为常数。由于异重流潜入受水流、泥沙及边界条件复杂，影响因素较多，现有的理论模型通常未考虑潜入点流速的分布，将动量修正系数的影响考虑为常数，虽然给出了异重流潜入条件中含沙量与潜入点弗劳德数的显性关系式，但对于动量修正系数的确定，仍缺乏理论依据。

由于不同水库自然条件千差万别，加之水库实际调度过程各不相同，目前尚未能对其水沙运动规律准确把握与模拟，使得对潜入点的认识还停留在经验性较强的阶段，现有的处理方法难以反映潜入点的流速分布状况和动量修正系数的变化情况。本文根据异重流潜入点流速分布特征，给出潜入点流速分布，得到了动量修正系数的理论结果。

1 潜入点垂线流速分布形式 1.1 潜入点位置

潜入点处漂浮物聚集，水面出现有回流现象，常成为判断潜入点位置的直观标志。异重流潜入位置随着来水来沙条件和水下地形的变化而变化。从平面上看，随着扩散角的增大，依次为主流半月形、单侧半月形、全断面半月形，这是由于异重流沿横向流速分布为二次曲线形，潜入主流中间部分流速大。当潜入处的断面过宽时，潜入后的异重流并不分布于整个库底，主要区别是水下异重流有双侧回流、单侧回流、无回流。异重流发生后，流动经过宽阔库段，异重流发生扩散，如图 1所示^[10]。综上所述，可以从直观上判断异重流潜入与否。随着水位高低、入库流量的大小及底部淤积情况等不断变化，潜入点位置不仅上、下游移动，而且潜入处的异重流也会左右横向摆动^[11]。

挟沙水流与水库内清水存在密度差，在一定的地形边界条件下，二者之间密度差引起的压力差导致浑水开始潜入库底，形成异重流，在潜入点附近的水面上，上游带来的大量漂浮物，在潜入点下游清水的倒流作用下聚集 (图 2、3)。

1.2 流速与含沙量分布

当水库发生异重流时，在潜入点附近，水流由普通浑水明流转化为异重流。由于水深沿程增加，其流速和含沙量沿垂线的分布形状将沿程发生变化^[12]，具体可描述如下：

1) 在离潜入点较远的上游，水深较小，流速较大，含沙量较大，流速和含沙量沿水深呈正常分布；

2) 到离潜入点不远处，水深增大，流速和含沙量分布呈不正常状态，最大流速位置向库底移动；

3) 在水深增大到一定程度，浑水开始潜入库底，此处为异重流潜入点，这里流速及含沙量沿垂线分布很不均匀，在水面处流速为0，含沙量也几乎为0，最大流速位置进一步向库底靠近；

4) 潜入点往下，异重流已经形成，异重流的流速和含沙量沿水深分布比较均匀，异重流之上形成横轴环流，含沙量的零点在水面以下。在潜入点处，有漂浮物聚集，这通常是判定异重流发生的一个直观标志。潜入点的水流泥沙条件可以作为判定异重流是否发生的条件。

从图 4可看出，明流转化为异重流时主流线流速垂线分布沿程的变化情况。此处尽管上层清水能被浑水带动，但是由于闸门开启，清水被带出库外，不能形成横轴环流，看不到表层清水的向上流动。

前人的研究认为，从异重流潜入点及其以下河段的流速垂线分布沿程变化看：在水库明流段，水库通过挟带水流与河床泥沙发生交换；在潜入点上游，由于壅水作用，水流表面流速由大变逐渐减小；至潜入点处，流速最大值出现在水下，水流开始分层，浑水潜入清水下层，在潜入点处的表面流速和底部流速均为零。根据水库测验资料和水槽试验资料，异重流的流速分布接近抛物线型，最大流速点发生在交界面以下，含沙量的分布则在交界面突然增大，然后比较均匀增加，直至底部。这也与钱宁^[13]的认识一致。

根据图 4中异重流流速分布形式，假设异重流流速分布接近抛物线形，服从表层及底层流速均为0近似的二次抛物线规律，即：

$ \frac{{{u}^{'}}}{{{V}_{\text{m}}}}=a{{\left( \frac{z}{{{h}_{\text{p}}}} \right)}^{2}}+b\left( \frac{z}{{{h}_{\text{p}}}} \right) $

(1)

式中：u′为异重流垂线上某一点z的流速, m/s; V_m为异重流的垂线上最大流速, m/s; h_p为异重流潜入点浑水水深，m。

对2001—2015年小浪底水库实测异重流潜入点进行了分析整理，根据流速分布图形，结合含沙量垂线分布特点，分成3类:1) 上层有薄层清水；2) 明显地分层，上下层流动方向相反；3) 其流速分布同明渠流流速分布类似。由于潜入点处水流变化剧烈，不易稳定测验，因此对实测数据进行了合理性分析。结果见图 5。

按照流速分布特点，当$ \frac{z}{{{h}_{\text{p}}}}=0, 1$时，$ \frac{{{u}^{'}}}{{{V}_{\text{m}}}}=0$时，$ \frac{z}{{{h}_{\text{p}}}}=\frac{1}{2}$，$ \frac{{{u}^{'}}}{{{V}_{\text{m}}}}=1$可得出式 (1) 的参数a=-4, b=4, 代入式 (1) 得：

$ \frac{{{u}^{'}}}{{{V}_{\text{m}}}}=4{{\left( \frac{z}{{{h}_{\text{p}}}} \right)}^{2}}+4\left( \frac{z}{{{h}_{\text{p}}}} \right) $

(2)

根据小浪底水库实测资料和物理模型试验数据拟合结果，可得a=-1.79, b=2.36, 代入式 (1) 得：

$ \frac{{{u}^{'}}}{{{V}_{\text{m}}}}=1.79{{\left( \frac{z}{{{h}_{\text{p}}}} \right)}^{2}}+2.36\left( \frac{z}{{{h}_{\text{p}}}} \right) $

(3)

可以看出，两个系数的理论值和实测值是有差别的，差别在于潜入点是一个由明流过渡到异重流的变化区域，受测验方法、阻力和水流、泥沙运动的复杂性的影响，导致理论值与模型试验、实测资料的拟合值的差别，但两者数值接近，正负及数量级相同。

2 潜入点判别条件 2.1 动量修正系数理论分析

建立的非恒定异重流运动的动量方程：

$ \begin{align} &{{J}_{\text{b}}}-\frac{{{f}_{\text{m}}}}{8}\frac{u_{\text{m}}^{2}}{{{\eta }_{\text{g}}}^{'}g{{h}_{\text{m}}}}-\frac{{{\tau }_{\text{c}}}}{{{h}_{\text{m}}}\left( {{k}_{\text{m}}}{{\gamma }_{\text{m}}}-{{\gamma }_{\text{c}}} \right)}-\frac{\partial {{h}_{\text{m}}}}{\partial x}= \\ &\ \ \ \ \ \frac{1}{{{\eta }_{\text{g}}}^{'}g}\left( \frac{\partial {{u}_{\text{m}}}}{\partial t}+{{\alpha }_{\text{m}}}{{u}_{\text{m}}}\frac{\partial {{u}_{\text{m}}}}{\partial x} \right) \\ \end{align} $

(4)

式中：u_m为异重流平均流速；τ_c为附加阻力；$ {{J}_{\text{b}}}=-\frac{\partial {{Z}_{b}}}{\partial x}, {{Z}_{\text{b}}}$为河底高程；α_m为因流速沿垂线不均匀分布引起的动量修正系数；f_m为异重流综合阻力系数；γ_m、γ_c为平均浑水容重、清水容重；h_m分别为浑水层的厚度；${{\eta }_{\text{g}}}^{'}=\frac{{{k}_{\text{m}}}{{\gamma }_{\text{m}}}-{{\gamma }_{\text{c}}}}{{{k}_{\text{m}}}{{\gamma }_{\text{m}}}}$为重力修正系数；g为重力加速度。式 (4) 中的动量修正系数α_m的计算公式为：

$ {{\alpha }_{\text{m}}}=\frac{\int_{0}^{{{h}_{\text{p}}}}{{{u}^{'}}{{u}^{'}}\text{d}z}}{u_{\text{m}}^{2}{{h}_{\text{m}}}} $

(5)

式 (4) 中的含沙量修正系数k_m的计算公式为：

$ {{k}_{\text{m}}}=\frac{\int_{0}^{{{h}_{\text{m}}}}{\left( \int_{z}^{{{h}_{\text{m}}}}{{{\gamma }_{\text{m}}}^{'}\text{d}z} \right)\text{d}z}}{0.5{{\gamma }_{\text{m}}}h_{\text{m}}^{2}} $

(6)

k_m的大小取决于含沙量S沿水深的变化趋势，目前可利用常见的含沙量沿水深分布公式近似计算。

对于恒定异重流，存在$ \partial {{u}_{\text{m}}}/\partial t=0$，因此式 (4) 中的水沙因子仅随x而变。在2维恒定流情况下，单宽流量 (q_m) 沿程不变，故应有：

$ \begin{align} &\frac{\text{d}}{\text{d}x}\left( {{q}_{\text{m}}} \right)=\frac{\text{d}}{\text{d}x}\left( {{h}_{\text{m}}}{{u}_{\text{m}}} \right)={{h}_{\text{m}}}\frac{\text{d}}{\text{d}x}\left( {{u}_{\text{m}}} \right)+{{u}_{\text{m}}}\frac{\text{d}}{\text{d}x}\left( {{h}_{\text{m}}} \right)=0 \\ &或\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \frac{\text{d}{{u}_{\text{m}}}}{\text{d}x}=-\frac{{{u}_{\text{m}}}}{{{h}_{\text{m}}}}\frac{\text{d}{{h}_{\text{m}}}}{\text{d}x} \\ \end{align} $

(7)

将式 (7) 代入式 (4) 化简可得：

$ \frac{\partial {{h}_{\text{m}}}}{\partial x}=\frac{{{J}_{\text{b}}}-\frac{{{f}_{\text{m}}}}{8}\frac{u_{\text{m}}^{2}}{{{\eta }_{\text{g}}}^{'}g{{h}_{\text{m}}}}-\frac{{{\tau }_{\text{c}}}}{{{h}_{\text{m}}}\left( {{k}_{\text{m}}}{{\gamma }_{\text{m}}}-{{\gamma }_{\text{c}}} \right){{\eta }_{\text{g}}}^{'}g}}{1-\frac{{{\alpha }_{\text{m}}}}{{{\eta }_{g}}^{'}g}{{u}_{\text{m}}}\frac{{{u}_{\text{m}}}}{{{h}_{\text{m}}}}} $

(8)

已有异重流潜入点形成条件的水槽试验观测结果表明^[3]：从明流过渡到异重流的清浑水交界面曲线上发现有一拐点k，其$ \text{d}{{h}_{m}}/\text{d}x\to-\infty $。同样采用该条件，则由式 (7) 可得：$ {{\alpha }_{\text{m}k}}u_{\text{m}k}^{2}/\left( {{\eta }_{g}}^{'}g{{h}_{\text{m}k}} \right)=1$。由于潜入点处在k点上游，潜入点水深h_p大于h_mk，则u_p显然比u_mk要小，即应有：

$ \frac{{{\alpha }_{\text{m}}}u_{\text{p}}^{2}}{{{\eta }_{g}}^{'}g{{h}_{\text{p}}}} < 1 $

(9)

利用式 (1) 分别积分计算可得：

$ {{u}_{\text{m}}}=\frac{\int_{0}^{{{h}_{\text{p}}}}{{{V}_{\text{m}}}\left( \frac{bz}{{{h}_{\text{p}}}}+\frac{a{{z}^{2}}}{h_{\text{p}}^{2}} \right)\text{d}z}}{{{h}_{\text{p}}}}={{V}_{\text{m}}}\left( \frac{b}{2}+\frac{a}{3} \right) $

(10)

$ \begin{align} &\int_{0}^{{{h}_{\text{p}}}}{{{u}^{'}}{{u}^{'}}\text{d}z}=V_{m}^{2}\int_{0}^{{{h}_{\text{p}}}}{\left( \frac{bz}{{{h}_{\text{p}}}}+\frac{a{{z}^{2}}}{h_{\text{p}}^{2}} \right)\text{d}z}= \\ &\ \ \ \ \ \ V_{\text{m}}^{2}{{h}_{\text{p}}}\left( \frac{{{b}^{2}}}{3}+\frac{{{a}^{2}}}{5}+\frac{ab}{2} \right) \\ \end{align} $

(11)

已知h_p=h_m, 将式 (10)、(11) 代入式 (5) 可得：

$ {{\alpha }_{\text{m}}}=\frac{\int_{0}^{{{h}_{\text{p}}}}{{{u}^{'}}{{u}^{'}}\text{d}z}}{u_{\text{m}}^{2}{{h}_{\text{m}}}}=\frac{\left( \frac{{{b}^{2}}}{3}+\frac{{{a}^{2}}}{5}+\frac{ab}{2} \right)}{{{\left( \frac{b}{2}+\frac{a}{3} \right)}^{2}}} $

(12)

当z=h_p时，理论值$ {{u}_{\text{m}}}=\frac{2}{3}{{V}_{\text{m}}}$, 代入式 (12) 可得：

$ {{\alpha }_{\text{m}}}=1.20 $

(13)

当z=h_p时，实测值u_m=0.34V_m, 代入式 (12) 可得：

$ {{\alpha }_{\text{m}}}=1.13 $

(14)

式 (12) 即为动量修正系数的理论表达式，从式 (13) 中可以看出，动量修正系数与潜入点水深成反比，水深越大，动量修正系数越小。对于潜入点处流速沿垂线分布形式为抛物线分布时，其理论值为式 (13)，对于小浪底水库，其计算值为式 (14)。

李书霞等^[7]认为流速沿垂线分布服从对数规律，则α_m=1+g/C²κ²，其中, C为谢才系数；κ为von Karman常数，一般取0.4。如流速垂线分布符合1/7幂函数规律，则α_m=1.016；如果服从表层及底层流速均为0近似的二次抛物线规律，则计算的动量修正系数α_m一般大于1.1。计算结果与其预测相符合。

2.2 潜入点判别条件推导

根据以往研究，在异重流潜入点的水流流态不仅与流速沿垂线不均匀分布引起的动量修正系数 (α_m) 有关，而且还与参数η_g′有关。而2.1节分析表明，η_g′取值与含沙量的大小、分布形式 (悬浮指标)，及压力不均匀分布系数k_m等都有关系。

利用浑水容重γ′、清水容重γ₀、泥沙干容重γ_s、含沙量S可得^[11]:

$ {{\eta }_{\text{g}}}'=\frac{{{\gamma }^{'}}-{{\gamma }_{0}}}{{{\gamma }^{'}}}=\frac{\left( {{\gamma }_{\text{s}}}-{{\gamma }_{0}} \right)S}{\left( {{\gamma }_{\text{s}}}-{{\gamma }_{0}} \right)S+{{\gamma }_{0}}{{\gamma }_{\text{s}}}} $

(15)

当γ_s=2 650 kg/m³时，式 (15) 可变为:

$ {{\eta }_{\text{g}}}'=\frac{1.65S}{1.65S+2\ 650}=1-\frac{1}{1+\frac{S}{1\ 606}} $

(16)

对式 (16) 分别进行一次求导和二次求导分别可得：

$ {{\left( {{\eta }_{\text{g}}}' \right)}^{'}}=\frac{\frac{1}{1\ 606}}{{{\left( 1+\frac{S}{1\ 606} \right)}^{2}}}>0 $

(17)

$ {{\left( {{\eta }_{\text{g}}}' \right)}^{''}}=-\frac{\frac{1}{1\ {{606}^{2}}}}{{{\left( 1+\frac{S}{1\ 606} \right)}^{3}}} < 0 $

(18)

从式 (17) 可知，重力修正系数η_g′随含沙量S单调递增，从式 (18) 可知，曲线是凹型曲线。这就表明重力修正系数η_g′随含沙量S的增大而增大，而且增大的幅度在逐渐减小。

根据含沙量S_v与η_g′关系可以得到：

$ {{\eta }_{\text{g}}}'=-1.747\ 4{{\left( {{S}_{\text{V}}} \right)}^{2}}+1.611\ 5{{S}_{\text{V}}}+0.001\ 1 $

(19)

已有水槽试验结果显示^[12]，异重流潜入点处的密度Froude数Fr_p′会随含沙量的增加而减小，大约含沙量为400 kg/m³时为止；且又在含沙量为400 kg/m³时发生转折，Fr_p′值会直线下降。以往研究通常根据水槽试验资料点绘Fr_p′与含沙量 (S) 或体积比含沙量 (S_V) 之间的相关关系。由于Fr_p′本身间接含有含沙量因子，显然这种处理方法是不合理的，需要进一步改进。

因此式 (8) 改写为：$ \frac{u_{\text{p}}^{2}}{g{{h}_{\text{p}}}}=f\left( {{\alpha }_{\text{m}}}, {{S}_{\text{V}}}, \frac{\omega }{\kappa {{u}_{*}}} \right)$。潜入点处流速沿垂线不均匀分布程度前文已进行了计算，取1.13。受目前实测资料的限制，还无法直接考虑悬浮指标对潜入点形成的影响。因此本次研究仅考虑体积比含沙量对异重流潜入点形成的影响，即：

$ Fr_{\text{p}}^{2}=u_{\text{p}}^{2}/\left( g{{h}_{\text{p}}} \right)=f\left( {{S}_{\text{V}}}/{{\alpha }_{\text{m}}} \right) $

(20)

式中：Fr_p为潜入点处的Froude数；f(S_V) 为体积比含沙量S_V的某一函数，需要由实测资料确定。

采用已有资料主要集中在多沙河流水库发生的异重流，如小浪底水库、官厅水库，测验时含沙量一般不会超过400 kg/m³，故大于400 kg/m³的水槽试验结果在本次分析中暂不考虑。点绘异重流潜入点处的水槽试验资料^[13-15]，如图 6所示，Fr_p与S_V之间的关系可用式 (21) 表示：

$ \begin{align} &Fr_{\text{p}}^{2}=u_{\text{p}}^{2}/\left( g{{h}_{\text{p}}} \right)=0.263{{\left( {{S}_{\text{V}}} \right)}^{0.76}}, \\ &或\ \ \ \ F{{r}_{\text{p}}}={{u}_{\text{p}}}/\sqrt{g{{h}_{\text{p}}}}=0.51{{\left( {{S}_{\text{V}}} \right)}^{0.38}} \\ \end{align} $

(21)

将式 (18) 与式 (20) 联立，可得潜入点处的密度Froude数为：

$ \begin{align} &F{{r}^{'}}_{\text{p}}^{2}=u_{\text{m}}^{2}/\left( {{\eta }_{\text{g}}}g{{h}_{\text{m}}} \right)=0.24{{\left( {{S}_{\text{V}}} \right)}^{-0.137}} \\ &或\ \ \ \ F{{r}_{\text{p}}}^{'}={{u}_{\text{m}}}/\sqrt{{{\eta }_{\text{g}}}g{{h}_{\text{m}}}}=0.49{{\left( {{S}_{\text{V}}} \right)}^{-0.069}} \\ \end{align} $

(22)

将q_m=u_ph_p代入式 (22)，则可得潜入点水深与单宽流量及体积比含沙量之间的关系表达式，即

$ {{h}_{\text{p}}}=0.738\ 6q_{\text{m}}^{2/3}/S_{\text{V}}^{0.764/3} $

(23)

式 (23) 表明:当来流单宽流量增加时，潜入点应下移；当来流含沙量增加时，潜入点应上移。在实际异重流潜入条件判别中，一般可以采用式 (23) 估算潜入点水深。

整理小浪底水库2006—2015年调水调沙期^[16]、三门峡水库1961—1962年异重流期间实测资料，代入式 (23) 进行估算，潜入点处计算与实测水深对比如图 7所示。从图 7中可以看出，计算与实测的潜入点水深较为接近。

原型观测数据中潜入点水深在3.8~21.0 m，流速在0.19~1.35 m/s，含沙量在4~169 kg/m³之间变化。

但正如范家骅^[13]指出的，由于现场情况的复杂性及测验仪器精度限制，原型观测资料大多不是真正潜入点处进行的测验，图中部分数据较为分散。因此可以认为式 (23) 可以用于小浪底水库异重流潜入点的判别条件。

3 结论

通过理论推导、实体模型试验与实测资料分析，对水库异重流潜入点流速分布形式公式进行了研究，得到以下结论：

1) 在对异重流潜入点流速垂线分布进行文献分析的基础上，提出其数学表达式近似接近截距等零的抛物线，抛物线的系数理论值分别为-4、4。将小浪底水库实测的异重流潜入点流速沿垂线分布资料进行率定，其系数分别为-1.79、2.76。两者数值接近，正负及数量级相同，其中的差别在于实际测验的误差、浑水阻力变化与水流、泥沙运动的复杂性。

2) 根据动量修正系数的定义，将潜入点流速垂线分布公式进行计算后，相应得到动量修正系数的理论值为1.2，小浪底水库实测及物理模型实测值为1.13，二者甚为接近。

3) 根据2) 的结果，收集2006—2015年小浪底库区潜入点、1961—1962年三门峡库区资料、小浪底库区模型试验资料，结合前人研究，对潜入点判别条件进行了修正，得到了新的潜入点修正公式，其计算结果与实测值接近。

建议进一步开展异重流潜入点的含沙量沿垂线分布研究，加强原型观测，提高异重流潜入点判别条件的理论性，增加对异重流潜入点的认识。