颗粒相质量、动量和脉动速度的二阶矩方程

颗粒相质量、动量和脉动速度的二阶矩方程

1关于s的扩散理论虽然悬浮沉淀颗粒浓度分布的理论多种多样，但每一种结果表明的形式都与扩散方程相似，只是一些不同的扩散系数。长期以来,悬沙浓度分布的研究一直以传统的泥沙扩散理论为主,泥沙扩散系数εs的确定是理论的关键。但是,泥沙扩散系数的确定往往依靠半经验的处理,假定εs等于或者正比于水流的紊动粘性系数vtf。根据已知泥沙浓度的垂线分布,利用传统扩散方程反推出εs往往大于vtf,这样的结果与离散质点的紊流扩散理论并不一致。虽然针对扩散系数的确定提出了种种修正方法,如Czernuszenko提出了漂移扩散系数的概念,但都未从根本上解决这个矛盾。而且,对于试验和实测中所观测到的泥沙浓度I型分布,传统的泥沙扩散理论亦难以解释。由于传统理论存在不足,一些学者着手寻求更为一般的理论解释悬沙运动机理。如Wang&Ni应用固液两相流的动理学理论,分析发现升力作用是I型分布形成因素之一。刘大有从两相流的一般扩散模型出发,发现传统扩散理论不能反映颗粒脉动强度梯度引起的扩散,指出在颗粒脉动强度变化较大区域,传统泥沙理论所基于的Fick扩散定律存在明显缺陷,并对εs大于vtf的试验结果以及浓度I型分布的形成做了定性解释。刘大有继而分析了泥沙悬浮运动的机理,认为河流或水平管道中颗粒的悬浮,除了重力和浮力以外,还可能受到其他的升力与颗粒脉动强度梯度的作用。本文从动力学理论的PDF方程出发,修正传统的泥沙扩散方程,并对悬沙浓度不同分布类型的形成机理及其影响因素作定量分析。2u3000颗粒相动与颗粒动脉内流动力学模型低浓度挟沙水流中,不考虑泥沙颗粒之间的碰撞,单颗粒运动的Lagrange方程为dRpidt=νpi(1)dνpidt=(ui-νpi)/τp+Fi+wi(2)式中,Rpi为颗粒位置,νpi为颗粒速度,ui为颗粒位置处的水流速度;τp为颗粒驰豫时间;Fi为阻力以外的其他相间力与重力之和,wi为颗粒布朗运动作用力,是一连续随机过程。将Lagrange方程转化成为Euler变量的方程,定义颗粒相平均的PDF(概率密度分布函数):〈f〉=〈δ(xi-Rpi)δ〈νi-νpi)〉(3)式中,〈·〉表示系综平均,δ(·)为Dirac的δ函数。颗粒相平均浓度〈C〉和平均速度〈νi〉:〈C〉=∫〈f〉d→ν,〈C〉〈νi〉=∫νi〈f〉d→ν(4)模化颗粒轨道的水流紊动与颗粒的作用〈u′iδ(xi-Rpi)δ〈νi-νpi)〉,对式(3)求时间导数,可导出平均PDF的输运方程:∂〈f〉∂t+νi∂〈f〉∂xi+∂∂νi〈f(〈ui〉-νiτp+Fi)〉=Gp〈u´iu´k〉∂2〈f〉∂xi∂νk+1τp{fp〈u´iu´k〉+Dδik/τp}∂2〈f〉∂xi∂νk(5)式中,D为颗粒布朗运动扩散系数,〈ui〉和u′i分别为水流平均与脉动速度;fp和Gp为描述颗粒与水流紊动作用的参数。不失一般性,采用:fp=1-exp(-ΤE/τp),Gp+fp=ΤE/τp(6)式中,TE为颗粒和水流紊动涡的作用时间。对方程(5)两边分别乘以1,→ν和→ν→ν,并对整个速度空间积分,可进而获得颗粒相连续方程、动量方程和脉动速度二阶矩方程:∂〈C〉∂t+∂〈C〉〈νi〉∂xi=0(7)∂〈νi〉∂t+〈νj〉∂〈νi〉∂xj+∂〈ν´jν´i〉∂xj=〈ui〉-〈νi〉τp+〈Fi〉-Dtpijτp∂ln〈C〉∂xj(8)∂〈ν´iν´j〉∂t+〈νk〉∂〈ν´iν´j〉∂xk+1〈C〉∂〈C〉〈ν´iν´jν´k〉∂xk=-Dtplk1τp(δilδjm+δjlδim∂〈vm〉∂xk+2τp(fp〈u´iν´j〉+Dδij/τp-〈ν´iν´j〉)(9)式中,ν′i=ν′-〈ν′i〉为颗粒相脉动速度;Dtpij=τp(〈ν′iν′j〉+Gp〈u′iu′j〉)为颗粒扩散张量,包括颗粒无规则脉动和水流紊动扩散两部分。3泥沙扩散方程明渠二维恒定均匀流中,∂/∂t=0、∂/∂x=∂/∂z=0,x为流动方向,z为侧向,y为垂向。颗粒相连续方程简化成:∂〈C〉〈νy〉∂y=0(10)对方程(10)积分得〈C〉〈νy〉=Jw=const.,由自由水面处颗粒相零通量条件Jw=0,可得〈νy〉=0,这里〈νy〉为Euler观点的颗粒相平均速度,并非单颗粒沉速;同理有液相垂向平均速度〈uy〉=0。颗粒相垂向动量和垂向脉动速度二阶矩〈ν´2y〉的方程分别为:∂〈ν´2y〉∂y=〈Fy〉-1τpDtpyy∂ln〈C〉∂y(11)∂〈C〉〈ν´3y〉∂y=2〈C〉τp(fp〈u´2y〉+D/τp-〈ν´2y〉)(12)方程(11)可作为研究悬沙垂线浓度分布的基本方程。为与传统的泥沙扩散方程相比较,引入挟沙水流中泥沙的沉速ωp,它由阻力、重力和浮力之间的平衡获得:ωp=τp(1-ρf-ρs)g(13)在垂向作用力中,除了颗粒所受有效重力(1-ρf/ρs)g以外,还包括Saffman升力和Magnus升力等,将这些升力表示成综合升力FL,则有:〈Fy〉=FL-ωp/τp(14)将式(13)和式(14)代入方程(11),得到二维明渠恒定均匀流的泥沙浓度方程:Dtpyy∂〈C〉/∂y=-ωp〈C〉+〈C〉τp(FL-∂〈ν´2y〉/∂y)(15)或者:Dtpyy∂〈C〉/∂y=-ωp〈C〉β-11(15*)β1={1-FL-∂〈ν´2y〉/∂y(1-ρf-ρs)g}-1(16)方程(15)左端反映了挟沙水流中泥沙颗粒的扩散,右一项为重力沉降,右二项为颗粒相脉动强度梯度和升力所产生的垂向合力的影响。β1表征了这种合力对有效重力作用的增强或削弱。如果忽略右二项,即β1=1,则方程(15)为传统的泥沙扩散方程:εs∂〈C〉/∂y=-ωp〈C〉(17)其中,εs表示传统扩散方程定义的扩散系数,以与Dtpyy区别。在后文中,不引起混淆,用C表示平均浓度〈C〉,用Dtp表示颗粒垂向扩散系数Dtpyy。4影响浓度分布的因素4.1细颗粒泥沙垂向及布朗运动Dtp浓度方程(15)中的颗粒垂向扩散系数Dtp和脉动速度二阶矩〈ν´2y〉需要封闭。对于细颗粒,局域平衡假设近似满足,或者假定ν´y服从正态分布,方程(12)中的扩散项可略去:〈ν´2y〉=fp〈u´2y〉+D/τp(18)这样,细颗粒泥沙的脉动将包括由于水流紊动引起的颗粒脉动和无规则的布朗运动,垂向扩散系数Dtp为:Dtp=ΤE〈u´2y〉+D(19)壁面附近的粘性底层中,TE〈u´2y〉趋于零,布朗运动扩散系数D为常数,这样,细颗粒泥沙的布朗运动将居主要地位,使得泥沙的浓度梯度维持有限值。粘性底层以外,TE〈u´2y〉远大于D,细颗粒泥沙将以紊动扩散为主,Dtp≈TE〈u´2y〉。因为细小颗粒的“穿越轨道效应”并不明显,颗粒与水流紊动作用时间尺度TE与载能涡的时间尺度Tt近似相等,则有:Dtp≈Tt〈u´2y〉=εm,εm为水流垂向紊动扩散系数,这与单相流体紊流扩散理论一致。4.2边壁周边区域实测资料表明,壁面附近以外的区域,〈ν´2y〉沿垂向接近均匀分布,或者呈现出下面大、水面小的分布形式,亦即∂〈ν´2y〉/∂y≤0。在壁面附近区域,由于水流垂向脉动强度向壁面递减,可能出现∂〈ν´2y〉/∂y>0。由于泥沙颗粒密度大,其沿流向的速度分量通常小于水流的速度分量,作用于颗粒上的升力将大于零,且在边壁附近显著。定义升力系数αL,u*为水流剪切流速,作用于细颗粒上的升力可表示为:FL=αLdpyu3*νfρs-ρfρs√νfu*y(20)定义无量纲垂向坐标η=y/H,径深比d0=dp/H,固液密度比s=ρs/ρf、沙粒雷诺数Rep*=u*dp/νf、二维水流能坡Je=u2*/gH,将式(20)代入式(16)有:β1={1-αL√d0Rep*Jeη3/2+sJe(s-1)u2*∂〈ν´2y〉∂η}-1(21)由上式可见,β1的取值将存在多种可能,在壁面附近以外区域,若∂〈ν´2y〉/∂η≤0沿垂向近似保持不变,升力和脉动强度梯度所产生的合力大于零。径深比d0、沙粒雷诺数Rep*、能坡Je、水沙密度比1/s的增加,都将使合力增大。在合力数值小于有效重力时,β1>1;等于有效重力时,|β1|=∞;大于有效重力时,β1<0。在壁面附近区域,由于∂〈ν´2y〉/∂y可能大于零,上述合力可能为正、也可能为负。当合力为负值时,有β1<1。水沙密度比1/s的增加,将促进合力负值和β1<1的出现。4.3参数l型的出现在传统泥沙运动理论中,见方程(17),用实测浓度分布资料反算出的泥沙扩散系数εs通常大于εm,亦即β1≥1。根据上文分析,比较方程(15)与(17)有:εs=-ωpC∂C/∂y=Dtpβ1,β=εsεm=εsDtp=β1(22)可见,用已知浓度分布反算泥沙扩散系数时,方程(17)忽略了升力和颗粒脉动强度梯度的影响,通常过高估计了泥沙颗粒的重力沉降作用,导致εs值较大,从而常有β1≥1。而实际上,由于β1<1在一定流动条件下也可能出现,β<1也是可能的。而且,由于εs取正值,方程(17)给出∂〈C〉/∂y<0,不能预测浓度分布I型的存在。方程(15)中,当有效重力小于升力和脉动强度梯度的合力时,β1取负值(此条件下的等式(22)没有意义,用传统扩散方程已不能定义εs),∂〈C〉/∂y>0,浓度I型分布出现。因此,方程(15)是传统扩散方程(17)的修正方程,β1为其修正系数。5u3000不同升力和颗粒相脉动强度梯度对浓度分布和传统扩散方程的影响粘性底层以外,泥沙扩散系数Dtp中布朗运动的贡献可以略去:Dtp≈Τt〈u´2y〉=εm≈νtf(23)对于低浓度、含细颗粒的明渠二维挟沙水流,液相水流紊动涡粒性系数νtf和垂向脉动强度〈u´2y〉可用下式计算:νtfu*Η=κ(1-η)[1η+πΠsin(πη)]-1(24)√〈u´2y〉=1.23u*exp(-0.67η)(y+=yu*/νf≥50)(25)式中κ为Karmann常数,取为0.41;Π为尾流强度。对于清水,Π=0～0.2,随水流雷诺数变化;对挟沙水流,Π与Richardson数有关,Coleman给出Π=0.2～0.9。通过Π的变化,式(24)一定程度上反映挟沙水流流速分布偏离清水流速分布、从而涡粘性系数不同的事实。而且,本文考察升力和颗粒相脉动强度梯度对浓度分布和传统扩散方程的影响,采用式(24)计算νtf对结果没有本质影响。对粒径dp=0.1～1.0mm的单颗粒,清水中的静水沉速ω0和挟沙水流中的沉速ωp:ω0=10νfdp⋅[(1+0.01(ρs/ρf-1)gd3p/ν2f)0.5-1](26)ωp/ω0=(1-C)5(27)用方程(15)、(18)、(21)、(23)～(27)可以计算分析不同因素对β1和用传统扩散方程预测浓度分布所带来的误差的影响。用BL和Bv分别表示升力和颗粒相脉动强度梯度与有效重力的比值、定义无量纲沉速ω*=ω0/u*、颗粒Stokes数St=τp/Tt,将方程(15)写成:∂C∂η=-ω*β1κC(1-C)m1+ηπΠsin(πη)η(1-η)(28)β1=(1-BL-Bv)-1(29)BL=αLη-3/2√d0Je√Rep*(30)Bv=1.2321-1/sJeexp(-1.34η)⋅{1.34+exp(-1/St)(S-2t∂St/∂η-1.34)}(31)式中的ω*/κ和ω*/β1κ即为通常的悬浮指标和修正悬浮指标。式(30)中的升力系数αL可用Saffman力估计,结果为αL～1.62/(3πk1.5)=0.7。实际上,除了Saffman力以外,颗粒还受其他升力的影响,计算中取αL=0.8。边界条件理论上可通过一定的颗粒-边壁作用模型建立,为简便,仍用参考点和参考浓度的方法,下面除特别说明以外,参考点ηa=0.05,参考浓度Ca=0.5%,尾流强度Π=0.2,计算域最低点y+=yu*/νf=50,最高点η=0.9。不同流动条件下的计算参数和典型位置处的β1的计算结果见表1,β1、方程(28)的解C(η)、β1=1时方程(28)的解C0(η)(即传统扩散方程的解)、以及ΔCrel=(C-C0)/C0的计算结果如图1～5所示。5.1不同剪切流速下的升力变化水深不变时,u*增大使得Je和Rep*增加,升力参数BL增大,而颗粒相脉动强度梯度参数Bv变化不大,故β1相应增大。从图1来看,β1随η的减小而增大,在η=0.1以上,不同剪切流速下的β1都不大;而在接近边壁时,β1迅速增大。这表明升力变化主要影响壁面附近域内的β1。在Je=0.8%时,η=0.008以下的β1出现了负值(图中没有绘出),对应的C(η)最大值点在床面以上,为I型浓度分布。ΔCrel随着剪切流速u*的增加而增加。Je=0.2%时,计算域底部和顶部的相对差别均不超过10%;而在Je=0.8%时,顶部ΔCrel为17%,底部达45%。4种剪切流速下,ΔCrel在参考点以上为正、参考点以下为负,表明C(η)比C0(η)分布均匀。5.2dp对rel活性的影响随着dp增大,Rep*、d0和St数增大。这样,尽管Bv稍微减少,但BL增大导致β1增大。由于ω*也随之增大。ΔCrel将随dp增大而迅速变化。从图2可见,在η=0.1以下,dp=0.1mm时的β1变化最小,其C(η)和C0(η)几乎重合,ΔCrel最大不超过6%;dp=0.3mm时的β1变化最大,最大值达3.4,ΔCrel在底部达40%。在η=0.1以上,dp=0.1mm时的ΔCrel〈3%,而dp=0.3mm时达22%。5.3为1.0.6crel的小密度比密度比s的变化,对BL没有直接影响。s变小使得主流区的Bv增加和β1增大。但在近壁区,当颗粒相脉动强度为正梯度时,Bv取负值,s的减小将导致β1的减小。从图3看,对于s=1.056和1.6的小密度比情形,随着η的减小,β1开始增大,随后减小,并甚至出现β1<1,这是由于Bv为负值,且绝对值较大所致。ΔCrel在小密度比s时都较小,s=1.056和1.6时的ΔCrel均在7%以内。而s=7.8时的ΔCrel在顶部达到26%,在底部达到30%。这是由于s小时,虽然在主流区β1较大,但ω*很小,导致浓度分布均匀。降低了C(η)和C0(η)的相对差别。5.4体大学生c、曲线随着水深H的增大,能坡Je和径深比d0减小,导致BL和Bv减小,β1随之减小并趋向于1,但ω*保持不变。这样,不同水深条件下的C0(η)相同,C(η)随着H的增大而趋近于C0(η),如图4所示(图中仅绘出了一个C0(η)曲线,其他情形与其重合)。图4表明,在η=0.1以上,H=0.5m和1.0m时的β1都接近于零;在同一η处,较小的H对应较大的β1。ΔCrel分布图上,H=0.05m时的顶部取值25%,底部取值也比同一η处的其他水深情形为大;而H=0.5m和1.0m时的顶部取值在2%以内,底部取值也不超过17%。5.5ca-ct对crel分布的影响参考点处的Ca从0.05%到3%变化。由于C的变化对Je、d0、Rep*、ω*没有影响,仅改变沉速ωp、Stokes数St和Bv,对β1和∂C/∂η构成一定影响。图5中,随着Ca的变化,β1仅在底部有一定变化,在η=0.05以上几乎没有变化。ΔCrel的分布在计算域内(y+>50和η<0.9)也几乎与浓度变化无关。这表明在低浓度下,浓度变化对ΔCrel近似没有影响。应该说明的是,随着浓度C的增加,式(24)中的尾流强度Π增大或者Karman常数κ减小,水流紊动涡粘性系数随着变化。而这里的分析中,并没有将这些变化考虑在内。6关于1和crel的垂线分布下面结合Wang&Qian和Einstein&Chien的细颗粒试验资料分析上述计算结果。各组次的试验及计算参数见表2。其中,尾流强度Π随着试验条件的变化有所调整。Wang&Qian的SF组次s小、Rep*大;Einstein&Chien的S11组次Je和Rep*都较大。图6为计算的β1和ΔCrel的垂线分布。由图可见,SF组次的β1在主流区较大,而在近壁区小于1;S11组次的β1在主流区并不明显,在近壁区较其它组次为大。这是由于密度比s的变化主要影响Bv,Je和Rep*的变化主要影响BL造成的。不同试验组次的ΔCrel相比较,由于SF1和SF2的ω*比其他组次小得多,在主流区尽管其β1较大,但其浓度分布均匀,因而ΔCrel最小;S11组次的ω*最大,浓度分布最不均匀,由于β1也较大,故其在整个计算域内的ΔCrel最大,顶部将近40%,底部在50%以上。图7为计算的C(η)和C0(η