摘要：近年来用于研究含沙水流中悬浮颗粒垂向浓度分布的理论已有很多。连续性假设虽然被证明在描述流体运动时非常成功，但却不足以描述含沙水流中的离散固体颗粒运动。随机模型能用于研究流动中单个颗粒的运动，但很难解释固体颗粒之间相互作用的机制。本文对各种传统理论进行了综合分析和比较，对已有的典型颗粒浓度分布的一般性解释进行了讨论，并据此提出了今后研究的重点。
关键词：水沙流 泥沙悬浮 连续介质理论
悬浮颗粒垂向浓度分布被认为是研究含沙水流中颗粒运动特性的主要指标。这项有意义的研究以Rouse经典理论的提出和随后Vanoni的实验研究为标志，并在此后取行了很大的进展[1，2]。许多学者提出了各种理论和公式。事实上含沙水流可以被看作一个两相流系统，其中的液相和固相遵守基本的守恒定律，各相之间由相间耦合作用而联系。固液两相流系统可用宏观或微观的方法进行描述，如连续理论或动理论[3～12]。本文在已有工作的基础上进一步对各种理论进行比较，并对已有泥沙浓度分布公式更广泛的概括形式进行讨论。
1 现有理论比较
关于悬浮颗粒浓度垂线分布规律的研究已有很多。在众多的理论和模型中，应用较广泛的有扩散理论、混合理论、两相流理论、随机理论、动理学以及相似理论。
扩散理论要求分散的颗粒对水流结构有较小的影响，这意味着该理论只适用于尺寸和比重较小的颗粒。扩散理论的运用一般基于质量守恒和均匀紊流。假定水流紊动扩散作用和颗粒重力作用达到平衡，则可得到一个简单的扩散方程[1]
ωC εs(dC/dy)=0
式中 C为距离床面任意高度y处的悬浮颗粒浓度；ω为颗粒沉降速度；εs为泥沙扩散系数，这里假设其等于清水紊流的动量交换系数ε。基于上述方程求解得到的Rouse公式(见表1)在应用上获得了巨大的成功，以致于许多后继的研究者认为在应用时仅仅需要对该理论进行简单的修正或改进即可。由于假设扩散理论在低浓度含沙水流及细颗粒条件下是有效的，所以修正主要集中在两个方面：一是对泥沙扩散系数εs的修正；二是对悬浮指标Z=ω/κu*或悬浮指标中的参数(如ω或κ)的修正。尽管一些研究者指出Rouse公式能通过对参数Z、εs或κ的简单线性修正而应用于更广的范围，但在一定范围以外这种做法是有缺陷的[3～14]。理由如下：(1)由实测数据得不到εs和ε沿整个垂线的线性关系(见图1)；(2)Zm(实测悬浮指标)比Zc(计算悬浮指标)小的设想并不总是正确，图2就是一个反例[15]；(3)认为κ值可变而对卡门常数进行修正的尝试是不可靠的，因为事实上Coleman已发现(如图3所示)κ是一个不变的常数[16]。实际上，任何在扩散理论框架内进行的修正无非是寻找泥沙扩散系数的表达式。由于扩散理论没有给出悬浮颗粒运动的动力学解释，为了能更深入地探讨悬浮颗粒的机理，许多学者致力于寻找更为普遍的理论。
| 图1 实测εs和ε的垂线分布[2]
Vertical distribution of εsand ε from the measured data |
| 
| 图2 Zm和Zc的关系[15]
Relation of Zm and Zc from the measured data |
|
表1 悬浮颗粒浓度垂线分布的代表性公式 |
Representative formulas for vertical sediment distribution |
|
序号 | 作者 | 公式 | 说明 |
|
1 | LaneKalinske | C/Ca=exp[-6κ(ω0/u*)(y/H-a/H)] | 为距床面距离a处的颗粒浓度 |
2 | AnanianGarbashian | C/Ca=exp[-(y/H-a/H)/AA] | AA=(0.0017v2/gH)[ρP-(1 KA)ρ/(ρP-ρ)] |
3 | Cai | C/Ca=[(C1 y/H)/(C1 a/H)]-C2 | C1=B/A，C2=κ(ω0/u*)/A |
4 | Velikanov | 
| Φ=H(1-y/H)ln(1 y/Δ)Δ为床面当量糙率 |
5 | KarimKennedy | 
| Γ=[1/(1 m)](1-y/H)(y/H)1-m |
6 | LaursenLin | C/Ca=exp[ω(1 1/m)/βumaxf(Iy/H-Ia/H)] | 
|
7 | Laursen | 
| |
8 | TanakaSugimoto | 
| |
9 | Barenblatt | 
| |
10 | Hunt | 
| 0.995＜B*＜1 |
11 | Rouse | 
| |
12 | Zagustin | 
| 
|
13 | ItakuraKish | 
| φ=f(ω/u*) |
|
混合理论假设流体和分散的悬浮固体颗粒可分别视为连续介质[17，18]。尽管流体和固体颗粒具有不同的密度、速度和其它特征，但混合流中液相和固相的体积浓度可以很容易确定。具体做法是对两相分别列出动量守恒方程，随后将之叠加得到固液相混合的总体方程，以消除两相相互影响的复杂项。混合理论架建于单一流体模型的基础之上，能克服扩散理论的部分缺陷并给出低浓度颗粒垂线分布的动力学机理。然而该理论在应用上存在的困难是如何确定混合流体的本构关系。McTigue采用Drew的方法，得到了一个对应于分层模式的运动方程[17]。然而，经典的紊流扩散方程仍被用作一个描述雷诺平均紊流的模型，这意味着该理论不可避免地存在类似于扩散理论中存在的
问题，即颗粒垂线分布的解决仍依赖于εs的确定。
两相流理论在处理这一复杂过程时强调流体和悬浮颗粒间的相互作用[19]。固液两相分别用一系列守恒方程进行描述，并以耦合项联系两相。方程的闭合依靠关于各相的本构关系模型、热力学状态方程或紊流附加方程和各相的描述。两相的耦合有不同的方法，如从紊流单颗粒运动方程修正得到，或者用其它的闭合模型[20，21]。然而，由于数学上的困难，很难获得方程的显式解。而且，在
目前的测量水平下各相间的差别和相互作用很难精确地测量。在研究粘性颗粒运动、非均匀颗粒运动或高浓度固液两相流时会遇到更多的困难。
随机理论的建立基于紊动流体和颗粒运动都具有随机特性的认识，它把流体中的悬浮颗粒运动看作是在液体表面到底壁之间范围内的随机徘徊[22，23]。随机理论的核心是确定颗粒在沿水深的垂线上任意高度和时刻向任意方向(如向上、向下或停留)运动的可能性，即概率密度。通常假定连续的随机跳跃符合马尔可夫过程。应用随机理论的困难在于确定颗粒从一个高度跳到另一高度的空间和时间步长[24]。迄今，随机理论的应用仍仅限于均匀紊流。大的步长不易满足均匀紊流的假设，而太小的步长又很难保证连续的两次跳跃是相互独立的，即是否符合马尔可夫过程。所以，用随机理论几乎不可能描述相间或颗粒间的相互作用。有趣的是，当取时间步长趋于零的极限时，由随机理论可以重新推导出传统的扩散方程[52]。因此，得到颗粒垂线浓度分布仍依靠εs的确定。
能量理论要求对流体和悬浮颗粒的平均流能和涡能作出精确的表述，而这正是固液两相流研究中最困难的[25～27]。Velikanov从能量平衡的观点出发，提出了所谓的“重力理论”，认为耗散的悬浮功(或悬浮能量)被用于支持紊流中重颗粒悬浮[28]。通过一些过分的简化，Velikanov推导出一个更加复杂的颗粒垂线分布公式(见表1)[28]。尽管近年来在能量理论基础上做了许多工作，对Velikanov的成果在某些方面也进行了改善和提高，仍然有许多细节需要进一步的探索[14，29]。一个著名的观点(如，Bakmeteff[25])认为含沙水流中用于支持颗粒悬浮的能量不是直接取自有效势能，而是从紊动能中获得。从这一观点出发，蔡树棠仿照周培源提出的涡量相似理论假设流体和颗粒紊动具有某种相似性，得到了与Rouse公式相似的颗粒垂线分布公式[30，31]。 | 
| 图3 κ与Ri(Richard数)的关系[16]
κ-Ri(Richardson Number)relation |
|
然而，存在的问题依然是能量事实上是怎样耗散的?悬浮能量在总的紊动能中占多大比例?这些问题还远没有解决。有关水流紊动能的公式可以在大部分流体力学教科书中找到[32，33]。然而目前关于维持悬浮颗粒的紊动方面的认识还很有限。由于对含沙紊流中能量耗散方式的认识不足，能量理论的实际应用受到很大限制。
上述理论从不同的角度对紊流中的颗粒悬浮给出了不同的理解。每种理论都有利有弊，并都得到一个相应的颗粒垂线分布公式。然而从这些存在的理论可以发现尽管采用不同的理论和数学处理，最后都得到或接近按扩散方程所得到的形式，主要的差别是扩散系数εs的表达有所不同[34]。这一仅从形式比较(而非实质探讨)得到的认识构成了对表1所示的由各种理论获得的不同公式进行综合和统一的基础[1，19，28，30，35～42]。
2 现有公式的综合 有趣的是，对表1中的每个公式，研究者都能找到一些实测数据来验证其正确性并显示其优于它式，而当选用其它的资料时则得到相符很差或者完全不符的结论。例如[43，44]，由图4可以看到，即使采用同样的悬浮指数，从不同公式得到的颗粒垂线分布也有很大的不同。因此，事实给我们提出了这样一些问题：为什么每个公式都有一些实测资料相符而与其它资料不符?是否可能或怎样综合现有公式以得到一个一般化的公式?对此，本文将在文献[43]和[45]的基础上进一步拓广讨论。 一般而言，颗粒扩散系数εs的表达式为 | 
| 图4 相同悬浮指数不同公式的比较(Z=0.5)
Comparison of different formulas with same suspension index(Z=0.5) |
|
| (1) |
这意味着获得εs的关键在于确定颗粒垂向运动的特征长度L(通常指所谓的颗粒平均自由程)，或者相应的特征时间tm，这正是迄今为止紊流研究中最困难的问题。这里假定颗粒的均方值 (或均方根 在靠近床底以外的大部分区域是可预测的。由于以前的研究表明，紊流波动强度 沿垂线几乎不变，所以如何表达L就成为(1)式中的要害所在[46，47]。在这一领域的研究已有很多[24，48，49]，却未能找到适用于不同颗粒属性的一般公式[24]。从以往的研究来看，L可用如下的简单关系来概括，即 |
| (2) |
| | |
这里A和B是两个待定系数；m1和m2为指数，两者都是颗粒属性和边界条件的函数；H是水沙流深度。一旦L被确定，则颗粒扩散系数也相应地被确定。对低浓度含沙水流，现已获得颗粒浓度垂线分布的统一公式[44]
| (3) |
在低浓度情况下，这里的指数a通常取零。公式(3)不含有卡门常数κ，但包含一个和悬浮指数Z=ω/κu*类似的参数。
目前的大多数公式都能从式(3)中推演出来，最简单的例子是当m1=m2=0，或L/H=A B时。然而只有在均匀紊流时这种条件才得以满足。这种情况下，如果设A B=1/6，那么可以直接得到Lane Kalinske公式(见表1)[35]。当用一个综合参数ＡA代替给定值1/6时，可推导出相似的Ananian Garbashian公式(见表1)[19]。表1列出了AA的表达式，其中U是垂线平均速度；ρP和ρ分别是颗粒和流体的密度；KA是一个经验常数，尽管认为L沿整个流深不变理由并不充分，但在一些河流上的实测结果满足上述公式[14]。
大多数含沙水流并非均匀紊流，因此公式(2)中的指数m1and m2几乎不可能同时为零。当考虑床面的推移质作用时，B值不能为零。在这种情况下有两种选择“m1=1和m2=0”或“m1=0和m2=1”。按照前者可推演出蔡树棠的相似理论公式(见表1)[30]。按照后者，可以得到一个类似于Velikanov由重力理论所得到的公式[28]，但要注意到公式(见表1)中的对数项变化很小。在研究细颗粒或粗颗粒时，上述近似处理并不能很好地描述颗粒的垂线分布，如图5所示。 一般来说，当研究较低浓度含沙水流中的悬浮颗粒时，颗粒在床面进行推移运动的影响可以忽略，也就是说B=0的假定是比较合理的。因此，以下的讨论就在此基础上进行，即仍有A，m1和m2三个参数需要确定。当m2=1，A=1/(1 m)以及m1=1-m时(这里m是一个指数，广泛应用于指数型速度分布公式，通常取值为1/6～1/7)，可重新获得与Laursen Lin(1954)公式一致的Karim Kennedy公式[42]。 | 
| 图5 Mihaylova[50]实测数据和Velikanov公式的比较[28]
Comparison between measured data and Velikanov formula |
|
那么，未知的参数能否进一步减少或简化?这种可能性似乎存在，象A=m1=1这种近似对大多数情况基本上都可以接受。在同样的这种处理的基础上，大部分著名公式的相似和差别都可以得到比较。下面将用参数n取代m2来讨论怎样从公式(3)推导出表1中所列出的各种公式。
当n=0时，可得到Laursen公式[40]；当n=0.5时，可推导出Tanaka Sugimoto公式[38]。
注意到2(1-η)(1-1-η)和(1-η)0.8几乎等价(这里η=y/H)，可推导出基于能量理论的Barenblatt公式[37]。引入新提出的参数B*，由于其定义为0.995＜B*＜1，使之同样符合上面的等价条件，Hunt推出了相似的公式[36]。
当取n=1时，可推出著名的Rouse公式[1]；取n=1.5时，从方程(3)推出的公式除了表面附近外在整个垂直深度范围非常按近于Zagustin公式[39]。
此外，当n=2时，注意到φ的平均值大约为1，可得到一个类似Itakura Kishi的公式[41]。
当n取正整数时，可得到一个级数解[12]
| (4) |
相似地，n取任何值时都可解得颗粒垂线浓度分布。对不同的流体和颗粒属性，参数n取值不同。
3 结语
上述关于现有泥沙颗粒垂线分布公式的一般性表述形式已经由作者从两相流基本理论出发重新获得。可以看出，过去此类研究多源于连续介质理论。然而，当撇开理论的具体内涵而直接比较其最终方程形式时却都接近得到与扩散方程类似的形式，并具体体现为泥沙扩散系数εs的差异方面。这也说明从各种理论出发得到的诸公式有可能用更一般的形式表达。
几乎所有的从不同理论得到的公式都仅仅是当参数n取某些特殊值时统一公式的特例。因此，参数n的变化规律成为解决颗粒垂线分布问题的关键。即使对同样的悬浮指数，从方程(3)也可得到不同的理论模型，但这些公式可能互有很大的差别，如图4所示。从方程(3)推得的浓度分布公式的总是随n变化，所以对不同公式修正悬浮指数或卡门常数是没有意义的，因为每个公式都相应于给定水沙条件下n取不同值时所对应的特定形式，所以修正Rouse公式没有通用性。目前，能用于研究n值变化的实测数据还非常有限[2，13，18，51]，参数n的确定还有赖于更多的系统观测资料。
减小字体
增大字体