如何计算沙粒阻力与沙波阻力

根据第二节中提到的沙粒阻力的概念，在水流作用于床面的全部剪切力τ0中，只有一部分直接作用于颗粒上，引发推移质的运动和床面形态的变化过程（沙波的形成），这就是所谓的“沙粒阻力”，它是水流与颗粒本身之间直接作用的结果。在床沙有运动、但没有沙波（动平整床面）的情况下，阻力完全由边界非滑移条件及颗粒绕流（包括水流直接作功使推移质颗粒运动）所产生，此时沙粒阻力等于全部阻力。床面上有沙波时，则由沙粒阻力和沙波阻力共同构成全部阻力。沙波阻力是沙波绕流所产生的阻力，又称形态阻力（或压差阻力，其起因是沙波迎、背水面所受动水压力不同），它所对应的是较大尺度的紊动（与沙波的尺度相仿）。沙波阻力引起的时均能量消耗对推移质运动没有贡献，但转化为紊动动能而支持了悬移质的运动。计算沙波阻力方法的关键，在于给出沙波阻力与水流条件的图解或函数关系，使阻力方程封闭。

（1）Einstein-Barbarossa方法。将阻力分成沙粒阻力和形状阻力的计算方法是首先由Einstein和Barbarossa（1952）提出的。其途径与Einstein分割边壁阻力和床面阻力的方法类似，是将水力半径（或过水断面面积）分为两部分。Einstein和Barbarossa认为水流中单位水体消耗的时均能量可按水力半径分割法（或能坡分割法）分成两部分，分别对应于沙粒阻力和沙波阻力。按照剪切应力的叠加可以写出

式中：分别是沙粒阻力和沙波阻力对应的剪切应力；分别是相应于这两者的水力半径。

这就是床面阻力计算中的Einstein水力半径分割法。显然

对于沙粒阻力，可由下列Manning-Strickler（Strickler，1923）公式或平均流速对数公式求得沙粒阻力对应的水力半径

式中：；χ是修正系数，对粗糙紊流，可取χ=1，式（3-83）需要试算求解。

对于沙波阻力，认为是与沿床面的输沙率有关的函数。于是根据野外实测资料，提出了一个经验性的沙波阻力与水流条件图解函数关系，即

如图3-26，式中，而Ψ是Einstein水流参数，定义式为

图3-26　沙波阻力和水流参数的关系图（Einstein和Barbarossa，1952）

可见Ψ是Shields数的倒数。当Ψ的范围在1.6～50之间时，水流处于低能态区，也就是说床面形态以沙垄为主。Ψ＜1.6时，水流处于过渡区至高能态区，沙垄将夷平，或出现逆行沙垄，所以沙波阻力较小，的值不大。由图3-26可见，这一方法的应用范围主要局限于低能态区。如果已知河流的J、D65、代表粒径D、U等水力要素和床沙资料，就可通过试算求得和由图解求得，从而计算得到平均水深H。

（2）Engelund方法。明渠水流阻力的实质就是时均机械能量转化为紊动能量而耗散，即总水头的损失。水头损失又可以分为沿程损失和局部损失两种。这一概念是Engelund方法的出发点（Engelund1966，1967），即床面阻力引起的水头损失hf可分解为河床表面的沙粒阻力引起的沿程损失和沙波形态引起的局部损失

对于二维流动的情况，Engelund认为能坡损失完全由床面阻力引起（即，不计边壁阻力引起的JW，令Jb=J）。设沙波的波长为L，在一个波长范围内的水力坡降为

由此可见，Engelund法是能坡的划分。式中J′是坡降中由沙粒阻力引起的部分，而J″是由局部突然放大损失引起。其中，用水力学中处理突然放大引起的水头损失的方法，可按图3-27的二维沙波概化图形求出沙波形状引起的局部损失为

其中α是系数。如果水流的流动和沙波形状都是二维的，并取单位河宽进行计算，即Rb=R，τb=τ0，则有A1=h—a/2，A2=h+a/2，其中A1是波峰处的过水断面面积，A2是是波谷处的过水断面面积。设单宽流量为q，利用水力学中的封闭管道突然扩大损失的Carnot公式，有

图3-27　Engelund的二维沙波概化图形

设h≫a，q/h=U，则上式可简化为

记沙波波长为L，则可得出与沙波损失对应的能坡为

因此根据沙粒阻力和沙波阻力对能坡的分割就可以写为

上式两边同乘以后，变化为无量纲形式。其中D是代表粒径，可取为D50，能坡分割的无量纲形式为

式（3-86）中的各项就是Engelund提出的无量纲剪切应力（或称水流强度），记为

综合无量纲剪切应力

沙粒阻力引起的无量纲剪切应力

沙波阻力引起的无量纲剪切应力

因而Engelund的阻力分割可以记为

Engelund提出了两条假设：

1）动力相似的河流其无量纲综合剪切应力Θ相等；

2）动力相似的河流沙波阻力损失占总阻力损失的比例相等。(www.xing528.com)

用下标“1”和“2”分别表示两条河流，则上述假设可记为：

1）Θ1=Θ2；

2）或者说

显然，对于动力相似的河流，据Engelund假设必有Θ1=Θ2和即，对于两条不同的河流，只要无量纲剪切应力Θ′相等，则无量纲剪切应力Θ必然也相等，Θ和Θ′两者是一一对应的。这表明存在沙波阻力与水流条件的函数关系Θ=f（Θ′）。这样只要河流动力相似，则不论是小河流或室内试验结果，还是大型河流的实测结果，表达成无量纲剪切应力后，均可点绘在同一条曲线上。而在Einstein的方法中，建立图解关系时只使用了河流的实测资料。

Engelund没有进一步给出计算式（3-86）和式（3-89）中的未知数J′、a的方法，也未求解函数关系Θ=f（Θ′），而是转向分析试验和实测资料，从中归纳得到Θ=f（Θ′）的经验关系。Engelund在分析资料时作了如下假定

即，Engelund仍沿用水力半径分割的方法，计算沙粒阻力所对应的剪切应力。其中沙粒阻力所对应的水力半径R′，系采用对数型断面平均流速公式试算求解，在式中Engelund取ks=2D65，并认为水流处于阻力平方区（在这种情况下河流才能够动力相似），χ可取为常数1.0，这样对数型断面平均流速公式可写为

式（3-90）与Einstein方法中所用的式（3-83）不同。由于在建立沙波阻力与水流条件的函数关系Θ=f（Θ′）时所用的是式（3-90），所以使用Engelund方法时也必须采用式（3-90），否则无法正确地使用他所提出的函数关系Θ=f（Θ′）。

求出R′后，用γR′J代替式（3-88）中的γh J′（在二维流动中认为h=R），就得到Θ′的值。Engelund由试验和实测资料中得到的Θ=f（Θ′）经验关系如图3-28。

Engelund和Hansen（1967）根据Guy等人（1966）的水槽试验资料，点绘出Θ=f（Θ′）的关系后得到如下沙垄区（低水流能态）、逆行沙垄区（高水流能态）拟合式，

图中高低水流能态区之间的间断，发生在Θ′约为0.55处。式（3-91）即是图3-28所示的两条直线，实际上不能包括所有试验点据。在随后的研究中，Engelund和Fredsoe（1982）又建议将式（3-91）中的曲线改为

图3-28　由试验和实测资料中得到的Θ=f（Θ′）关系（Engelund 1967）

Brownlie（1983）将式（3-91）的曲线扩展到Θ′＞1的高水流能态区范围，得到

乐培九等（1992）在验证长江阻力时，将上述关系概化为

（3）Engelund方法的进一步完善。Engelund的Θ=f（Θ′）经验关系是依据D50=0.19，0.27，0.39，0.93mm的一系列动床水槽试验资料而建立的，因此适用范围也不可超出这一粒径范围之外。此外，Engelund的Θ=f（Θ′）经验关系，没有明确给出低水流能态到高水流能态的过渡区内关系。而在我国的长江中游、黄河下游，水流条件处于低水流能态到高水流能态的过渡区，是常见的现象。为了在更广泛的范围内应用Engelund方法，王士强（1990）、王士强和White（1993）把资料扩展到了0.018mm至28.65mm的水槽试验及天然河渠阻力资料，并着重研究了过渡区的阻力关系。试验中Rb/D50的范围为10～120000，分析中引入了一个新的参数D*=D50[g（γs/γ—1）/ν2]1/3，以表达粒径的影响，并以D*为参数对这些试验结果和实测资料进行了适线拟合，如图3-29，从图中资料拟合的曲线方程为

图3-29　考虑粒径影响和过渡区情况的Θ=f（Θ′）经验关系

图3-30　水流能态、床面形态判别图解

低水流能态区（曲线3′，4′，6′）：

式中：x=lg（Θ′/0.04）；lgK1=0.513—0.123lgD*—0.141（lgD*）2；lgK2=0.56—0.0647lgD*—0.2183（lgD*）2；lgK3=0.017—0.0347lgD*—0.2728（lgD*）2。

式中：ν为水的运动粘滞系数。

当D*＞80以后，Θ与Θ′的关系变化很小，可按常数D*=80计算。

高水流能态区（曲线1至6）

式中：A=1.4/lg（/0.04）；G=1+4.874exp（—0.79D*）；=0.68+0.32exp（—0.1D*）。

由图3-30可见，对应于每一个D50，关系曲线Θ=f（Θ′）在过渡区都是双值的，仍然要判断高低能态。为此王士强给出了用沙粒Froude数判断的临界线

对于已知Rb、U和D50的任意流动，如果其沙粒Froude数Fr*=U/小于式（3-98）给出的临界值，则流动处于低水流能态，床面形态为沙纹或沙垄。反之，流动处于高水流能态，床面形态为平整或逆行沙波，如图3-30。

（4）Kennedy-Alam-Lovera方法。Lovera和Kennedy（1969）仿照管道流动阻力和定边界明槽流动阻力的研究思路，点绘了水槽试验及野外实测的动平整床面沙粒阻力，得到了Darcy-Weisbach阻力系数f与Reynolds数的关系，如图3-31。可见，由于推动颗粒运动需要付出额外的机械能，动平整床面沙粒阻力系数f显示出随Reynolds数增大而增大的趋势，而不是象定床粘贴相同粒径的沙粒粗糙所形成的阻力那样，在Reynolds数较大时f成为一个常数。

图3-31　Lovera和Kennedy（1969）的动平整床面沙粒阻力系数

沙波阻力可以表达为形状阻力系数。为了求得沙波阻力的图解函数关系，Alam和Kennedy（1969）采用能坡分割方法，假定J=J′+J″。在量纲分析的基础上，得到形状阻力系数是沙粒Froude数和相对糙率Rb/D50的函数（如图3-32），即

图3-32　Alam和Kennedy（1969）的动床沙波阻力系数

由前面对不同单元阻力合成方法的分析可知，有沙波情况下的床面总阻力系数f是由f′和f″直接相加而得，f=f′+f″。对于已知横断面尺寸、中值粒径D50和比降J的冲积河流，可以根据图3-31和图3-32试算求解出其水位～流量关系。