顆粒相質量、動量和脈動速度的二階矩方程

顆粒相質量、動量和脈動速度的二階矩方程

1關于s的擴散理論雖然懸浮沉淀顆粒濃度分布的理論多種多樣，但每一種結果表明的形式都與擴散方程相似，只是一些不同的擴散系數(shù)。長期以來,懸沙濃度分布的研究一直以傳統(tǒng)的泥沙擴散理論為主,泥沙擴散系數(shù)εs的確定是理論的關鍵。但是,泥沙擴散系數(shù)的確定往往依靠半經驗的處理,假定εs等于或者正比于水流的紊動粘性系數(shù)vtf。根據(jù)已知泥沙濃度的垂線分布,利用傳統(tǒng)擴散方程反推出εs往往大于vtf,這樣的結果與離散質點的紊流擴散理論并不一致。雖然針對擴散系數(shù)的確定提出了種種修正方法,如Czernuszenko提出了漂移擴散系數(shù)的概念,但都未從根本上解決這個矛盾。而且,對于試驗和實測中所觀測到的泥沙濃度I型分布,傳統(tǒng)的泥沙擴散理論亦難以解釋。由于傳統(tǒng)理論存在不足,一些學者著手尋求更為一般的理論解釋懸沙運動機理。如Wang&Ni應用固液兩相流的動理學理論,分析發(fā)現(xiàn)升力作用是I型分布形成因素之一。劉大有從兩相流的一般擴散模型出發(fā),發(fā)現(xiàn)傳統(tǒng)擴散理論不能反映顆粒脈動強度梯度引起的擴散,指出在顆粒脈動強度變化較大區(qū)域,傳統(tǒng)泥沙理論所基于的Fick擴散定律存在明顯缺陷,并對εs大于vtf的試驗結果以及濃度I型分布的形成做了定性解釋。劉大有繼而分析了泥沙懸浮運動的機理,認為河流或水平管道中顆粒的懸浮,除了重力和浮力以外,還可能受到其他的升力與顆粒脈動強度梯度的作用。本文從動力學理論的PDF方程出發(fā),修正傳統(tǒng)的泥沙擴散方程,并對懸沙濃度不同分布類型的形成機理及其影響因素作定量分析。2u3000顆粒相動與顆粒動脈內流動力學模型低濃度挾沙水流中,不考慮泥沙顆粒之間的碰撞,單顆粒運動的Lagrange方程為dRpidt=νpi(1)dνpidt=(ui-νpi)/τp+Fi+wi(2)式中,Rpi為顆粒位置,νpi為顆粒速度,ui為顆粒位置處的水流速度;τp為顆粒馳豫時間;Fi為阻力以外的其他相間力與重力之和,wi為顆粒布朗運動作用力,是一連續(xù)隨機過程。將Lagrange方程轉化成為Euler變量的方程,定義顆粒相平均的PDF(概率密度分布函數(shù)):?f?=?δ(xi-Rpi)δ?νi-νpi)?(3)式中,〈·〉表示系綜平均,δ(·)為Dirac的δ函數(shù)。顆粒相平均濃度〈C〉和平均速度〈νi〉:?C?=∫?f?d→ν,?C??νi?=∫νi?f?d→ν(4)模化顆粒軌道的水流紊動與顆粒的作用〈u′iδ(xi-Rpi)δ〈νi-νpi)〉,對式(3)求時間導數(shù),可導出平均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為顆粒布朗運動擴散系數(shù),〈ui〉和u′i分別為水流平均與脈動速度;fp和Gp為描述顆粒與水流紊動作用的參數(shù)。不失一般性,采用:fp=1-exp(-ΤE/τp),Gp+fp=ΤE/τp(6)式中,TE為顆粒和水流紊動渦的作用時間。對方程(5)兩邊分別乘以1,→ν和→ν→ν,并對整個速度空間積分,可進而獲得顆粒相連續(xù)方程、動量方程和脈動速度二階矩方程:??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〉)為顆粒擴散張量,包括顆粒無規(guī)則脈動和水流紊動擴散兩部分。3泥沙擴散方程明渠二維恒定均勻流中,?/?t=0、?/?x=?/?z=0,x為流動方向,z為側向,y為垂向。顆粒相連續(xù)方程簡化成:??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)可作為研究懸沙垂線濃度分布的基本方程。為與傳統(tǒng)的泥沙擴散方程相比較,引入挾沙水流中泥沙的沉速ω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)為傳統(tǒng)的泥沙擴散方程:εs??C?/?y=-ωp?C?(17)其中,εs表示傳統(tǒng)擴散方程定義的擴散系數(shù),以與Dtpyy區(qū)別。在后文中,不引起混淆,用C表示平均濃度〈C〉,用Dtp表示顆粒垂向擴散系數(shù)Dtpyy。4影響濃度分布的因素4.1細顆粒泥沙垂向及布朗運動Dtp濃度方程(15)中的顆粒垂向擴散系數(shù)Dtp和脈動速度二階矩〈ν′2y〉需要封閉。對于細顆粒,局域平衡假設近似滿足,或者假定ν′y服從正態(tài)分布,方程(12)中的擴散項可略去:?ν′2y?=fp?u′2y?+D/τp(18)這樣,細顆粒泥沙的脈動將包括由于水流紊動引起的顆粒脈動和無規(guī)則的布朗運動,垂向擴散系數(shù)Dtp為:Dtp=ΤE?u′2y?+D(19)壁面附近的粘性底層中,TE〈u′2y〉趨于零,布朗運動擴散系數(shù)D為常數(shù),這樣,細顆粒泥沙的布朗運動將居主要地位,使得泥沙的濃度梯度維持有限值。粘性底層以外,TE〈u′2y〉遠大于D,細顆粒泥沙將以紊動擴散為主,Dtp≈TE〈u′2y〉。因為細小顆粒的“穿越軌道效應”并不明顯,顆粒與水流紊動作用時間尺度TE與載能渦的時間尺度Tt近似相等,則有:Dtp≈Tt〈u′2y〉=εm,εm為水流垂向紊動擴散系數(shù),這與單相流體紊流擴散理論一致。4.2邊壁周邊區(qū)域實測資料表明,壁面附近以外的區(qū)域,〈ν′2y〉沿垂向接近均勻分布,或者呈現(xiàn)出下面大、水面小的分布形式,亦即?〈ν′2y〉/?y≤0。在壁面附近區(qū)域,由于水流垂向脈動強度向壁面遞減,可能出現(xiàn)?〈ν′2y〉/?y>0。由于泥沙顆粒密度大,其沿流向的速度分量通常小于水流的速度分量,作用于顆粒上的升力將大于零,且在邊壁附近顯著。定義升力系數(shù)αL,u*為水流剪切流速,作用于細顆粒上的升力可表示為:FL=αLdpyu3*νfρs-ρfρs√νfu*y(20)定義無量綱垂向坐標η=y/H,徑深比d0=dp/H,固液密度比s=ρs/ρf、沙粒雷諾數(shù)Rep*=u*dp/νf、二維水流能坡Je=u2*/gH,將式(20)代入式(16)有:β1={1-αL√d0Rep*Jeη3/2+sJe(s-1)u2*??ν′2y??η}-1(21)由上式可見,β1的取值將存在多種可能,在壁面附近以外區(qū)域,若?〈ν′2y〉/?η≤0沿垂向近似保持不變,升力和脈動強度梯度所產生的合力大于零。徑深比d0、沙粒雷諾數(shù)Rep*、能坡Je、水沙密度比1/s的增加,都將使合力增大。在合力數(shù)值小于有效重力時,β1>1;等于有效重力時,|β1|=∞;大于有效重力時,β1<0。在壁面附近區(qū)域,由于?〈ν′2y〉/?y可能大于零,上述合力可能為正、也可能為負。當合力為負值時,有β1<1。水沙密度比1/s的增加,將促進合力負值和β1<1的出現(xiàn)。4.3參數(shù)l型的出現(xiàn)在傳統(tǒng)泥沙運動理論中,見方程(17),用實測濃度分布資料反算出的泥沙擴散系數(shù)εs通常大于εm,亦即β1≥1。根據(jù)上文分析,比較方程(15)與(17)有:εs=-ωpC?C/?y=Dtpβ1,β=εsεm=εsDtp=β1(22)可見,用已知濃度分布反算泥沙擴散系數(shù)時,方程(17)忽略了升力和顆粒脈動強度梯度的影響,通常過高估計了泥沙顆粒的重力沉降作用,導致εs值較大,從而常有β1≥1。而實際上,由于β1<1在一定流動條件下也可能出現(xiàn),β<1也是可能的。而且,由于εs取正值,方程(17)給出?〈C〉/?y<0,不能預測濃度分布I型的存在。方程(15)中,當有效重力小于升力和脈動強度梯度的合力時,β1取負值(此條件下的等式(22)沒有意義,用傳統(tǒng)擴散方程已不能定義εs),?〈C〉/?y>0,濃度I型分布出現(xiàn)。因此,方程(15)是傳統(tǒng)擴散方程(17)的修正方程,β1為其修正系數(shù)。5u3000不同升力和顆粒相脈動強度梯度對濃度分布和傳統(tǒng)擴散方程的影響粘性底層以外,泥沙擴散系數(shù)Dtp中布朗運動的貢獻可以略去:Dtp≈Τt?u′2y?=εm≈νtf(23)對于低濃度、含細顆粒的明渠二維挾沙水流,液相水流紊動渦粒性系數(shù)νtf和垂向脈動強度〈u′2y〉可用下式計算:νtfu*Η=κ(1-η)[1η+πΠsin(πη)]-1(24)√?u′2y?=1.23u*exp(-0.67η)(y+=yu*/νf≥50)(25)式中κ為Karmann常數(shù),取為0.41;Π為尾流強度。對于清水,Π=0～0.2,隨水流雷諾數(shù)變化;對挾沙水流,Π與Richardson數(shù)有關,Coleman給出Π=0.2～0.9。通過Π的變化,式(24)一定程度上反映挾沙水流流速分布偏離清水流速分布、從而渦粘性系數(shù)不同的事實。而且,本文考察升力和顆粒相脈動強度梯度對濃度分布和傳統(tǒng)擴散方程的影響,采用式(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和用傳統(tǒng)擴散方程預測濃度分布所帶來的誤差的影響。用BL和Bv分別表示升力和顆粒相脈動強度梯度與有效重力的比值、定義無量綱沉速ω*=ω0/u*、顆粒Stokes數(shù)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)中的升力系數(shù)α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。不同流動條件下的計算參數(shù)和典型位置處的β1的計算結果見表1,β1、方程(28)的解C(η)、β1=1時方程(28)的解C0(η)(即傳統(tǒng)擴散方程的解)、以及ΔCrel=(C-C0)/C0的計算結果如圖1～5所示。5.1不同剪切流速下的升力變化水深不變時,u*增大使得Je和Rep*增加,升力參數(shù)BL增大,而顆粒相脈動強度梯度參數(shù)Bv變化不大,故β1相應增大。從圖1來看,β1隨η的減小而增大,在η=0.1以上,不同剪切流速下的β1都不大;而在接近邊壁時,β1迅速增大。這表明升力變化主要影響壁面附近域內的β1。在Je=0.8%時,η=0.008以下的β1出現(xiàn)了負值(圖中沒有繪出),對應的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數(shù)增大。這樣,盡管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變小使得主流區(qū)的Bv增加和β1增大。但在近壁區(qū),當顆粒相脈動強度為正梯度時,Bv取負值,s的減小將導致β1的減小。從圖3看,對于s=1.056和1.6的小密度比情形,隨著η的減小,β1開始增大,隨后減小,并甚至出現(xiàn)β1<1,這是由于Bv為負值,且絕對值較大所致。ΔCrel在小密度比s時都較小,s=1.056和1.6時的ΔCrel均在7%以內。而s=7.8時的ΔCrel在頂部達到26%,在底部達到30%。這是由于s小時,雖然在主流區(qū)β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數(shù)St和Bv,對β1和?C/?η構成一定影響。圖5中,隨著Ca的變化,β1僅在底部有一定變化,在η=0.05以上幾乎沒有變化。ΔCrel的分布在計算域內(y+>50和η<0.9)也幾乎與濃度變化無關。這表明在低濃度下,濃度變化對ΔCrel近似沒有影響。應該說明的是,隨著濃度C的增加,式(24)中的尾流強度Π增大或者Karman常數(shù)κ減小,水流紊動渦粘性系數(shù)隨著變化。而這里的分析中,并沒有將這些變化考慮在內。6關于1和crel的垂線分布下面結合Wang&Qian和Einstein&Chien的細顆粒試驗資料分析上述計算結果。各組次的試驗及計算參數(shù)見表2。其中,尾流強度Π隨著試驗條件的變化有所調整。Wang&Qian的SF組次s小、Rep*大;Einstein&Chien的S11組次Je和Rep*都較大。圖6為計算的β1和ΔCrel的垂線分布。由圖可見,SF組次的β1在主流區(qū)較大,而在近壁區(qū)小于1;S11組次的β1在主流區(qū)并不明顯,在近壁區(qū)較其它組次為大。這是由于密度比s的變化主要影響B(tài)v,Je和Rep*的變化主要影響B(tài)L造成的。不同試驗組次的ΔCrel相比較,由于SF1和SF2的ω*比其他組次小得多,在主流區(qū)盡管其β1較大,但其濃度分布均勻,因而ΔCrel最小;S11組次的ω*最大,濃度分布最不均勻,由于β1也較大,故其在整個計算域內的ΔCrel最大,頂部將近40%,底部在50%以上。圖7為計算的C(η)和C0(η