扩散系数

悬浮颗粒在明渠剪切紊流中扩散系数计算公式研究 李 洪，许唯临，李克锋，李 嘉 (四川大学高速水力学国家重点实验室) 摘 要：本文利用各向同性均匀紊流理论推导了明渠剪切紊流中心区固-液两相速度差的计算公式，结果表明固液两相速度差依赖于颗粒密度、尺寸以及流体的平均速度和水力半径。进而研究了悬浮颗粒在紊流中的扩散系数，给出了紊流中的悬浮颗粒扩散系数的计算公式。该公式充分考虑了悬浮颗粒粒径、密度等颗粒特征参数以及流体粘性系数、平均流速和水力半径等水流参数对悬浮颗粒在紊流中的影响，能够比较全面地描述悬浮颗粒在紊流中的扩散情况。为了验证固-液两相流速度差计算公式，作者进行了水槽实验，实验结果与计算结果吻合良好。 关键词：固液两相流；扩散系数；速度差 基金项目：教育部科学技术重点项目资助 作者简介：李洪(1963-)，男，河南人，水力学及河流动力学博士生。 水利工程中的泥沙输移、环境工程中含固体颗粒的污水排放等都是典型的固液两相流问题，这些问题的求解在工程实际中具有重要的意义。悬浮颗粒在紊流场中的运动受控于紊流场、颗粒惯性和颗粒的自由沉降速度(交叉轨迹效应)，这些因素对颗粒的扩散都有影响，从而导致了颗粒扩散问题的复杂性。目前在求解河流中的泥沙输移、含固体颗粒的污水排放等固液两相流时，往往转换为对扩散方程的求解。颗粒态污染物由于吸附体的密度与水流密度不同，很难完全跟随水流一起运动。在目前的研究中，通常假定悬浮颗粒与水质点同步，因此不论是哪种颗粒，扩散系数都取同样的数值，并认为物质扩散系数与单纯流体的动量扩散系数相等，即Dp=Df，或是乘以某个系数以示区别。但是，实际上Dp=Df的假定仅适用于尺寸较小或密度接近于流体的颗粒。对于较大和较重的颗粒，这一假定已不适用。 为了推求颗粒的扩散系数，国内外都进行了大量的研究。对于定常均匀各向同性紊流场，Taylor[1]给出了时变的颗粒扩散系数的计算公式： Dpij(τ)= (1) 式中：Δuj(t)为第j方向t时刻颗粒相对于其平均值的速度差，〈〉表示系综平均。 由于无法考虑颗粒与流体点轨迹的不一致性，没有一种方法对具有有限尺寸和动量的真实颗粒的扩散给出了足够的描述。大部分作者都是假定Fick定律有效，并用Taylor公式计算颗粒的扩散系数。 Csanady[2]检验了重颗粒与流体点扩散的差别。假定颗粒的惯性足够小，可跟随流体质点的脉动，而其重量足够大，以致可产生一个自由沉速。Csanady分别用V0t/l1和Vgt/l两个数描述了垂直方向由于涡消失和交叉轨迹效应而引起u0p11(t)的变化。通过选择一个类似于Eulerian空间和Lagrangian流体点相关的up11(t)的函数形式，Csanady获得了颗粒纵向扩散系数Dp11(∞)的计算公式： Dp(11)(∞)=Df(∞){(Vg/V0)2+1/r2}-1/2 (2) 式中：Vg为定常漂移速度，V0为流体紊动强度，r=V0TL/l1，TL为 Lagrangian积分时间尺度，l1为垂直积分长度尺度。同样，他还获得了平行于Vg方向的颗粒扩散系数： Dp22(∞)=Df(∞){4(Vg/V0)2+1/r2}-1/2 (3) 可见Vg以倒数的比例形式减小颗粒的扩散系数。从上式可看出长时间的颗粒扩散系数与惯性无关，且大于相应的流体点的扩散系数，通常也大于基于Stokes阻力计算的颗粒扩散系数。由于颗粒雷诺数控制了附加于颗粒时间尺度的高阶项的作用，因而导致了颗粒扩散系数的增加。 Tchen[3]通过假定颗粒始终包含于一流体涡中，避免了求解颗粒运动方程的困难，这一结果被广泛引用。后来，由于“交叉轨迹效应”等直观概念的提出，澄清了造成颗粒流体扩散系数之间差别的两个机理：惯性的影响产生了颗粒扩散增大的趋势，颗粒从一与流体紧密相关的区域进入另一区域(交叉轨迹效应)产生了阻碍颗粒扩散的趋势。 Lee和Durst[4]，对颗粒在紊动管流中的运动作了详尽的研究，提出了颗粒频率响应扩散方法。Lee和Durst[5]认为，颗粒对于特征长度小于其直径(le≤dp)的涡旋将不作响应，它的运动完全由特征长度大于颗粒直径(le≥dp)的涡旋所支配。当le≥dp时，颗粒运动完全由紊动涡旋的扩散所支配，Dp=Df；当le≤dp时，颗粒运动独立于紊动涡旋的扩散，完全由周围流场平均运动的准层流粘性相互作用控制，Dp=0;当le≈dp时，颗粒运动既不完全由流体涡旋的扩散所控制，也不完全受周围流体的平均运动所支配。 本文拟通过对悬浮颗粒在紊流中跟随性的研究，给出一简单计算颗粒扩散系数的表达式。 1 本文计算方法 上述研究为避免和求解时的困难，要么假定紊流场为各向同性均匀紊流场，要么忽略Basset力的影响[8]。这在很大程度上限制了研究结果的推广应用。即使被广泛引用的Casandy公式，在推导时亦作了颗粒对流体具有良好跟随性的假定。关于Basset力对颗粒扩散的影响，Ahmadit[6]的研究表明：当扩散时间很长时，Basset力的作用不再重要。Gitterman[7]的研究表明：即使颗粒的密度与流体接近，颗粒的扩散系数亦会受到颗粒微尺寸的影响。 颗粒在流体中的跟随性取决于颗粒尺寸、密度以及水流条件等诸多因素，因此，颗粒扩散系数中应充分反映这些因素的影响。故计算颗粒扩散系数采用Hinze[8]的公式比较合理，Hinze认为速度的幅值比与颗粒的扩散特性直接有关，并提出了如下计算公式： Dp/Df= (4) 式中：Ef(ω)为角频率为ω的Lagrangian能谱，η是颗粒速度与流体质点速度的幅值比。但上式中包含了流体的Lagrangian能谱，而对于任意流动，要准确地给出流体能谱函数Ef(ω)的计算式是困难的。 当流动中的特征涡频率ωe可以确定时，即可计算出对应的ηe，则颗粒的扩散系数可按下式计算。 Dp=η2eDf (5) 一般要给出流动的特征涡频率ωe也是困难的。但ηe可由固液两相平均速度的比来计算，颗粒与流体速度的幅值比ηe可写为： ηe=1-r/uf (6) 式中f和r=f-p分别为流体速度和颗粒相对速度(两相速度差)的平均值，up为颗粒速度。考虑颗粒在主流向(纵向)的运动，并略Basset力项，则颗粒运动方程式可写为： (ρp+kmρf)dup/dt=(ρf+kmρf)duf/dt+18ρfνCD0/d2p(uf-fp) (7) 则颗粒相对速度满足下列方程： dur/dt=ρp-ρf/ρp+kmρfduf/dt-18ρfνCD0/d2p(ρp+kmρf)ur (8) 式中：CD0是颗粒阻力系数；dp为颗粒直径；ρf为流体密度；ν是流体运动粘性系数，经验常数km依赖于颗粒加速度的模数Ac，并有下列形式： km=1.05-0.066/(Ac+0.12) (9) 式中： Ac=|uf-up|2/dp(duf/dt-dup/dt) (10) duf/dt实际上是湍流中涡的加速度，与涡的尺度有关。设Tλ为颗粒以相对速度ur穿越尺度为λ的涡时所需的时间，将dur/dt理解为颗粒穿越该旋涡时旋涡对该颗粒相对速度变化的影响，则近似有Tλ=λ/ur及dur/dt≈ur/Tλ[9]，式(8)可改写成： dur/dt=Bτ/τ+Tλduf/dt (11) 式中：τ=d2p(ρp+kmρf)/18νρfCD0为颗粒的弛豫时间，B=ρp-ρf/ρp+kmρf。因dur/dt≈ur/Tλ=u2r/λ，于是有： u2r=Bτ/τ+Tλduf/dtλ (12) 由均匀各向同性紊流理论可知：当涡的尺度大于粘性耗散尺度λ0时，duf/dt=(ε2/λ)1/3，而涡的速度uλ=(ελ)1/3，其中ε为湍流能量耗散率。代入式(12)得： u2r=Bτ/τ+Tλu2λ (13) 上式即为两相速差的基本方程。 从式(13)中可看出ur依赖于涡的速度，进一步假定流动尺度不太小，颗粒尺寸不太大，这样对于大部分尺度λ＞λ0的涡，有τ≤Tλ。在此条件下，以Tλ代替式(13)中的τ+Tλ，于有是： u2r≈Bτ/Tλu2λ=Bτu2λur/λ (14) ur≈Bτu2λ/λ (15) 将τ的表达式代入有： ur≈(ρp-ρf/ρf)d2p/18νCD0u2λ/λ (16) 对于明渠流动，颗粒相对雷诺数一般都处于1＜Rep<1000的范围，取CD0=0.85Rep0.35，并把Rep=urdp/ν代入CD0，从式(16)得： ur≈0.144(ρp-ρf/ρf)1/1.35dp1.22/ν0.48(u2λ/λ)1/1.35 (17) 上式中，与涡尺度有关的部分是(u2λ/λ)1/1.35，因此速度差的统计平均值可表示为： r≈0.144(ρp-ρf/ρf)1/1.35dp1.22/ν0.48 (18)   其中表示(u2λ/λ)1/1.35的统计平均值。   由于采样过程的随机性，的表达式难于由直接的数学推导得出，应从实验中求得。忽略颗粒存在对流体的影响，则组成的一系列uλ和λ应仅与流动的水力半径R、流体的运动粘滞系数ν和流体流动的平均速度f有关，因此可表达成以下形式： ～ (19)   其中的x、y、z应通过实验和量纲分析求得。本文通过对9种颗粒在4种水流条件下得出的多组实验数据的平滑曲线可得到∝f1.62，再通过量纲分析可得： ∝f1.62/ν0.14R0.60 (20) 将上式代入式(18)可得： r≈k(ρp-ρf/ρf)1/1.35dp1.22f1.62/ν0.62R0.60 (21) 其中系数k由实验求得：k=0.035 这样，就得了明渠剪切湍流中心区颗粒与流体相对速度的近似表达式： =0.035(ρp-ρf/ρf)1/1.35dp1.22f1.62/ν0.62R0.60 (22) 路展民等[9]在直径为D垂直玻璃管道中的实验研究得出的两相速度差的表达式为： =0.0323(ρp-ρf/ρf)1/1.35dp1.24f1.58/ν0.58D0.66 (23) 上两式略有差别，这是由于所考虑的流动不同所致。 于是，由式(5)、(6)和(22)可得： Dp/Df=[1-0.035(ρp-ρf/ρf)1/1.35dp1.22f0.62/ν0.62R0.60]2 (24) 式中：ρp、dp、ν为颗粒和流体的特征参数，f和R也很容易测出，因此应用上式来计算颗粒的扩散系数比较方便。 对于明渠流，Df的计算公式为： Df=0.067Hu* (25) 将上式代入式(24)，可得： Dp=0.067Hu*[1-0.035(ρp-ρf/ρf)1/1.35dp1.22f1.62/ν0.62R0.60 (26) 上式即为计算悬浮颗粒在明渠中扩散系数的公式。 2 计算结果 为了检验公式(22)，作者计算了9种不同颗粒在不同水流条件下的相对速度，并与实验进行了对比，结果见表1。从表上可看出，固—液两相速度差强烈地依赖于流体的平均速度和颗粒的粒径和密度。利用本文公式(22)计算出的速度差与实测值符合很好，该公式充分反映了颗粒密度、尺寸和流体平均速度、粘性系数、水力半径对两相速度差的影响。 表1 不同水流条件下颗粒和流体的两相速度差(cm/s) f=16.0 31.5 42.5 71.5 颗粒密度/(g/cm3) 颗粒粒径/(mm) r测 r算 r测 r算 r测 r算 r测 r算 1.44 0.80～0.90 0.40 0.46 1.4 1.5 2.3 2.4 5.5 5.7 1.44 0.60～0.71 0.37 0.33 1.2 1.1 2.0 1.7 4.5 4.1 1.44 0.50～0.60 0.31 0.27 1.1 0.90 1.6 1.4 3.7 3.3 1.25 0.80～0.90 0.25 0.30 1.2 1.0 1.5 1.6 3.6 3.7 1.25 0.60～0.71 0.28 0.22 0.90 0.70 1.2 1.1 3.0 2.7 1.25 0.40～0.50 0.18 0.14 0.60 0.50 1.0 0.70 2.1 1.7 1.05 0.80～0.80 0.13 0.09 0.30 0.30 0.40 0.50 1.2 1.1 1.05 0.60～0.71 0.10 0.07 0.30 0.20 0.30 0.30 1.0 0.80 1.05 0.40～0.50 0.06 0.04 0.20 0.10 0.40 0.20 0.70 0.50 为了比较颗粒的扩散系数与流体扩散系数的差异，作者对前述9种颗粒情况进行了计算，结果见表2。 表2 不同水流条件下各种颗粒的扩散系数(×10-2cm/s) f/(m/s) 颗粒密度/(g/cm3) 颗粒粒径/mm 0.160 0.315 0.425 0.715 1.44 0.80～0.90 0.940 0.916 0.890 0.847 1.44 0.60～0.71 0.957 0.942 0.922 0.889 1.44 0.50～0.60 0.965 0.955 0.935 0.909 1.25 0.80～0.90 0.960 0.948 0.929 0.899 1.25 0.60～0.71 0.971 0.961 0.948 0.926 1.25 0.40～0.50 0.982 0.974 0.968 0.953 1.05 0.80～0.90 0.988 0.981 0.981 0.968 1.05 0.60～0.71 0.991 0.987 0.987 0.976 1.05 0.40～0.50 0.994 0.994 0.994 0.986 可见，在同一种流动中，对于同一种密度的颗粒，随着粒径的增大，颗粒扩散系数的减小增大；对于相同尺寸的颗粒，密度越接近流体，颗粒扩散系数的越接近流体扩散系数。 3 结束语 本文利用各向同性均匀紊流理论推导了明渠剪切紊流中心区固 液两相速度差的计算公式。结果表明：固液两相速度差依赖于颗粒密度、尺寸以及流体的平均速度和水力半径。本文将悬浮颗粒在紊流中的扩散系数与悬浮颗粒的跟随性和水流体条件结合起来，给出了紊流中的悬浮颗粒扩散系数的计算公式。该公式充分考虑了悬浮颗粒粒径、密度等颗粒特征参数以及流体粘性系数、平均流速和水力半径等水流参数对悬浮颗粒在紊流中的影响，能够比较全面地描述悬浮颗粒在紊流中的扩散情况。限于悬浮颗粒在流体中受力分析的复杂性，本文仅对球形颗粒在紊流中的扩散系数进行了研究，对于形状不规则的颗粒，更为深入的研究还需考虑颗粒的形状系数等多种因素的影响。