固体颗粒的群体沉降速度分析

1、固体颗粒的群体沉降速度分析郑邦民1，夏军强2(1.武汉大学 河流系，湖北 武汉430072; 2.清华大学 水利系，北京100084)摘要：从流体力学原理出发，数值模拟非均匀沙随机分布对流场的影响，推导出固体颗粒群体沉速的理论解。该公式不仅量纲和谐，浓度变化不超过极限浓度值，能反映含沙量与非均匀沙级配变化对群体沉速的影响，而且可避免其它公式量纲不和谐，计算中出现负值或降得过快的缺点。采用黄河实测资料对该公式进行了验证，计算结果与实测资料基本符合。关键词：固体颗粒; 群体沉速; 干扰流核；极限浓度1 引言泥沙在静止的清水中等速下沉时的速度，称为泥沙的沉降速度。在多沙河流的浑水中，泥沙颗粒的沉降特

2、性比清水中与低含沙水流中复杂。此时泥沙颗粒下沉相互干扰，部分颗粒或全部颗粒成群下沉，其下沉速度称为群体沉速1,2。群体颗粒沉降特性的研究具有十分重要的意义，它在多沙河流的河床演变分析和泥沙数学模型计算中广泛应用。单个颗粒的沉速与群体沉降可以相差10倍，故50年前有人说泥沙运动严格地讲只有一个半理论。为此应进一步分析颗粒群体沉降规律，使其在实际应用中不致有太大的误差。    本文在研究流体力学粘性流中圆球绕流规律的基础上，得出固体颗粒群体沉速的理论解，它可反映泥沙浓度与组成对群体沉速的影响。然后将该公式与现有的群体沉速公式进行比较，并用黄河实测资料进行验证。 2 理

3、论前提Navier_Stokes方程是流体力学的基本控制方程，它是求解流体力学诸多问题中普遍应用的方程。对不可压缩粘性流体，在有势外力作用下，可得Helmholtz涡量方程(1)上式中为流速矢量：为哈密顿算子(Hamilton Operator)； 为流体的运动粘滞系数;t为时间。一般情况下，三维流函数为向量，它与流速有如下关系。而流速与涡量，亦呈旋度关系，即。为了便于数值计算，它可写作一般曲线坐标系的张量形式：。其中。式中ui为逆变分量，j为协变导数，为协变基向量，它不一定是正交基，也不一定为单位基。对正交曲线坐标，则有其中uk为单位正交(局部)基上的物理分量；Hk为Lami系数或标量因子，

4、它反映微元弧长ddi与坐标微元di之间的比，即dsi=H(i)di。根据上述关系，我们可以将涡量方程写作一般曲线坐标形式或正交曲线坐标形式，以便于数值计算。它可以用来计算形体绕流等外部流动。对于二维流或在柱坐标、球坐标下的球对称，轴对称流动，式(2)可以简化。例如，在球坐标下有H1=1、H2=R、H3=Rsin，可得ds1=dR、ds2=d、ds3=Rsind。上式中R、为球坐标系下的三个坐标线。因轴对称时，且物理量只在R、方向上有变化，故有3=。同时可得R、坐标线上的速度分量 (3)这对小雷诺数下的圆球绕流，上述沉降分析是合适的。恒定流情况有惯性项均可忽略。对于外部绕流，流函数是无源场，则有

5、，因此可得。此时流函数与涡量的关系方程。如果是轴对称流动，则涡量只有3=为标量，流函数亦只有3=也为标量 ，而有+2=0              (4)    即在小雷诺数时，对轴对称的圆球绕流，解Navier-Stokes方程，可变为解流函数满足的重调和方程220。3 单个球形颗粒在粘性流中匀速沉降解单个球形细颗粒在粘性流中匀速沉降的速度0，可从流体力学分析得到。单个球形颗粒在粘性流中绕流时，其Stokes流函数为=1/4Vsin2

6、（3/R3R+2R2）。如将球坐标原点放在球心，利用(3)式，可得圆球绕流时R、坐标线上的流速分量分别为(5-1)(5-2)通过对作用于球面上的压力积分，可求得圆球所受阻力为3dV，其中V为球与流体的相对速度，当球体均匀沉降时，有效重力(s-)d3/6与阻力相平衡。其中球体半径为，直径为d。一个球体直径为d所占的距离为d+l，N个均匀颗粒占的距离当N(d+l)。一个球体体积为d3/6，N个d3/6，所占空间为N3(d+l)3，体积比浓度(6)当ld时，则Sv0；当l0，均匀沙排列均匀，得Sv=0.5236。此为极限浓度Svm的下界，随机紧密填充可达Sv=0.5612，如果为非均匀沙随机排列，该

7、值还可以再取高一些。如Svm0.65，但只要达到这种情况，流体将很难在颗粒间流动，因此，此下极限浓度值也是可用的，随着l/d的改变，浓度值变化如表1所示。表1 浓度Sv随l/d变化Table 1 Concentration Sv change with the variable l/d l/d 1005010752.710.50.30.2Sv0.5×10-50.4×10-54×10-40.0010.0240.010.06550.1550.23830.303不论如何，只要我们随机地给出粒径大小d与它所在位置，我们可以求得其它函数及阻力值。因为对于Stoke

8、s解，可以按奇异子线性叠加而得，可数值求解。而对于过渡区及紊流区，则非理论可解，而由实验决定。随着浓度v的增加，颗粒沉速有由过渡区趋向滞流区，紊流区趋向过渡区的趋势，因此重点放在理论分析滞流区沉降是合适的。 4 现有的群体沉速公式目前，对单颗粒泥沙在静水中的沉降规律己基本掌握，但对群体颗粒的沉降规律还有待于深入研究。前人对颗粒群体沉速公式的研究，可大致划分为两类：一是粗颗粒均匀沙的沉速，二是含较多细颗粒的非均匀沙沉速1。(1)Batchelor(1972)认为球体在低含沙水体中沉降时，颗粒间及颗粒与周围水体的相互影响，其沉速与其在无限清水中沉速的差异，是平均值不为0的随机变量3。他从统计理论出

9、发，最后推导出低含沙量情况下群体沉速的理论公式s/0=16.55Sv         (7)上式中当Sv0.05时，计算结果能与实验值基本符合；当Sv较大则偏差大。(2)Richardson和Zaki采用量纲分析与试验结果，建立如下群体沉速公式4 s/0=(1-Sv)m        (8)上式中指数m与沙粒雷诺数(Red=0d/)有关。夏震寰和汪岗对细沙取m=7时，上式与试验资料符合较好5。(3)王尚毅认为式(8)中当Sv=1时s=0，这种计算

10、结果不对6。因此将上式修改为s/0=(1-Sv)m    (9)上式中m=2.5；与泥沙特性有关，对塘沽淤泥可取=5.0。(4)钱意颖等人认为群体沉速的减小主要由于浑水的容重与粘度变化所致，得出了适用于层流区的群体沉速公式7(10)上式中、s、m分别为清水、泥沙及浑水的容重。(5)万兆惠等人认为细的单颗粒泥沙在清水中下沉时有(s-)d3/6=3d00。当为浑水时，上式仍成立，不过应以m代替0,m代替，s/(1-Sv)代替0。如浑水粘度采用日本森氏公式m/0=1+3Sv/(1-Sv/0.52)，代入上式可得群体沉速公式8s/0=(1-Sv)2/1+3Sv/(1-Sv

11、/0.52)            (11)(6)沙玉清认为在层流区，主要是浑水的粘度影响泥沙沉速，因此可得如下群体沉速公9(12)上式中d50取mm，且对d50在0.010mm附近的非均匀沙适用。(7)费祥俊认为用非均匀沙的中值粒径或平均粒径作为代表粒径，按均匀沙方法计算非均匀沙的平均沉速，将会导致较大的误差。因此应按各粒径组泥沙所占的比例，加权平均后得到非均匀沙的平均沉速公式10(13)式中Pi为第di粒径组泥沙所占的比例。浑水粘度m与含沙量大小和极限含沙量有关11。(8)张红武在

12、沙玉清公式基础上，考虑到沉降过程中一部分清水将依附沙粒同时下沉，结合试验结果，经推导得出如下群体沉速公式(14)上式中d50同样取mm。但该式适用范围比沙玉清公式大，近些年多用之于黄河泥沙数学模型计算。经数值计算我们发现，用式(14)计算群体沉速，必须使。5 本文的研究结果我们认为群体沉降公式在理论上要尽量合理，尽可能地有严格的两相流体力学的依据，量纲上要和谐，同时计算结果要与实测资料基本符合，才可用于实际计算。对于本文提出的群体沉速公式，作以下分析与论证。5.1 颗粒表面流速的分析从泥沙颗粒在浑水中受力情况进行分析：细颗粒泥沙沉降时阻力符合Stokes公式的单个颗粒沉降规律，为此多个颗粒的阻

13、力解是可以叠加的，只要是散粒体。我们可以在计算机上，做出随机变化的有限多个(10121015个)泥沙颗粒，粒径为0.100.01mm不均匀随机分布的泥沙颗粒受流体力的作用，从而得到由于泥沙下沉对周围流场的影响，这一影响并非简单地打一个(1-Sv)的折扣，而是对于周围流场的干扰，改变流函数、流线疏密形状的结果。    由流速u公式(5)中，可以看出：当R=a处，u=0，当R时，=90°时，u=V，这说明颗粒对流场有干扰，颗粒扰动形成流核，远处R>10时，u=V=0，流速等于沉速，球体匀速沉降，而当l不太大时，对u有一定影响。例如：当l=2,R+l=3

14、时，即s/0=0.74；当l=时，R+l=2时，即s/0=0.577；当l=0.1,R+l=1.1a时，则s/0=0.131。这些都说明浓度Sv的影响实质是对流场的影响。流体被干扰的流核，使其在一定柱状范围，要带动一定量的流体运动，其相对运动速度(沉速)降低了。这一结果反映于(15)式中。5.2浓度对群体沉降的影响  浓度对颗粒沉降影响的研究，最早是1906年A.Einstein从Brown运动得到一阶近似的理论结果，即m/=1+2.5Sv12。1972年Batchelor等人得到浓度影响的二阶近似理论结果(球体散颗粒)，相对粘性r=m/0=1+2.5Sv+7.6S2v。本文认为低浓

15、度 时的粘度改正应小些，高浓度时的粘度改正大些，非线性二阶式优于指数关系。图1给出了各家相对粘度公式的对比结果，可以看出Batchelor的二阶式居中。     考虑到群体沉降的极限浓度Svm及非均匀沙的颗粒组成的影响，本文给出的泥沙颗粒的群体沉降公式有如下形式(15)  图1相对粘度r与体积比浓度Sv关系(牛顿体)Fig.1Relationship between relative viscosity r and volumetric concentration Sv上式中为一修正系数，与混合沙的非均匀程度有关，对均匀沙，可取=1。

16、当Sv较小时，s对0改正不大，这是合理的，且s/0值大些。在中等浓度Sv下，则s/0偏差大些，在高浓度Sv时，各家差别更大，如(12)式、(14)式很快降为零。细粒泥沙的0本来就很小，取s/0=0在实用上是不方便的。上式中的极限浓度Svm，可采用方红卫确定的黄河干支流各站的结果13，也可用费祥俊提出的公式计算11。图2各家公式s/0的计算结果Fig.2Calculated s/0 from different formulas5.3各家公式比较与分析    图2给出了式(15)、式(10)、式(12)、式(14)在不同浓度下均匀沙(d=0.030mm)的群体 沉速

17、的变化规律。由图可知，沙玉清公式(式12)与张红武公式(式14)的计算结果较为接近。而本文提出的计算公式(式15)与钱意颖等人(式10)的计算结果相差不多，但式(15)考虑了极限浓度的影响，这比较符合实际情况。若采用式(15)计算非均匀沙的群体沉速，取d25=0.018mm、d50=0.030mm、d75=0.042mm，计算结果见图2。在相同的中值粒径和浓度下，泥沙组成偏细，采用式(15)计算的群体沉速可比式(10)小。这是因为细颗粒增多，导致悬浮液体粘性增大，从而使群体沉速降低。因此式(15)也可反映悬沙组成对群体沉速的影响。从公式形式上看，沙玉清公式(12)及类似公式(14)存在量纲不和

18、谐的问题。式(12)中，如要量纲合理，应当是，其中为无量纲系数；V为流速，可用沉速0代替，g为重力加速度。沙玉清利用明兹及赵乃熊试验结果分析了d50=0.10mm10mm时均匀颗粒的群体沉速与非均匀沙对比。结果表明当d50小时s/0改正多些；d50大时s/0变化小些。 当d50=0.1mm到10mm变化时，系数在301.33间变动。这样大的变化是难以选用的。公式(14)亦有此问题。关于粒径d的影响，一般已考虑在0中。对于均匀沙而言，滞流区有0d2，紊流区有，过渡区则介于两者之间变化。如浓度Sv之改变，使s与0在同一区内，则无d之影响。如Sv之改变使s由紊流区变为过渡区或滞流区，则s/0中出现d

19、的因素。但也不能认为其因子就是。早年 沙玉清提出的公式(12)及后人在“沙”的基础上提出的公式(14)都没有给出理论证明。式(7)(14)反映出均匀沙在同流区s/0的变化，无d之影响。不同区则s/0=dr，指数r=1.50之间，有了沙式(12)与张式(14)，不难用曲线拟合方式，给定系数与指数。因此，不一定需要在s/0式中采用如这种形式，因为不一定是唯一的，只是近似的经验表述(拟合)。它们都可能表示，随着d50之变细，s/0减小这一定性的变化规律。天然河道或水库中的泥沙，多为非均匀沙，其级配影响于之中，d25或d75越小，则s/0降得越多，反之亦然。5.4本文公式与实测资料的比较本文采用黄河支

20、流水文站60年代的群体沉速实测资料，对式(15)的计算精度作进一步的分析。选用的水文站分别为渭河的咸阳站、华县站，泾河的张家山站，北洛河的NFDA4头站。四站共有882组有效的实测数据，各站的悬移质含沙量大小及组成变化，如表2所示。表2各站的悬沙特征变化Table 2Characteristics changes of suspended load at different hydrological stations水文站有效数据(组)含沙量S/kg5m-3中值粒径d50/mm华县站 4110.14040.0030.0590.701.51咸阳站2250.14670.0030.0370.791.

21、50张家山站590.15180.0070.0730.821.18头站 1870.08810.0040.0520.761.18    在此采用费祥俊提出的方法，按各粒径组泥沙所占的比例，加权平均后得到非均匀沙的平均沉速公式。根据实测资料的率定，取式(15)中的参数=1.76。图3给出了实测与计算群体沉速的对比结果。由图可知，式(15)的计算结果与实测值基本符合。由图可知，在泥沙组成较细，含沙量较大时，群体沉速较小，计算值与实测值符合较好。当泥沙组成较粗时，计算值与实测值存在一定的误差。6结论本文从流体力学的基本原理出发，得出非均匀沙的群体沉速公式，该公式具有如下特点：反映了浓度(或含沙量)的影响；浓度不超过极限浓度值；引入反映非均匀沙组成对群体沉速的影响，它具有几何平均之意，既反映粒径的不均匀性，又具 有平均的代表性；计算结果能与黄河实测资料基本一致。图3实测与计算群体沉速的对比结果Fig.3Comparison between observed and c