第2第2查露等:三維顆粒群沉降的格子Boltzmann角度還是從基礎(chǔ)研究角度研究流體與懸浮顆粒之間相互作用都非常重要因此吸引了眾多研究］．在過去的幾十年中研究顆粒沉降的過發(fā)現(xiàn)顆粒尾渦對顆粒的運動產(chǎn)生重要的影響作用，其中最重要的是由2顆粒間相互作用而發(fā)生的拖曳－碰撞－翻滾(KT)運動．這種相互作用最早由Fotesetal．］試驗發(fā)現(xiàn)．在被發(fā)現(xiàn)之后，KT運動被眾多的研究復(fù)現(xiàn)包括實驗研究］，直接數(shù)值模擬(NS)6］以及格子oltmann()模擬］．限體積的顆粒因此研究實際懸浮顆粒與流體的相互作用是非常復(fù)雜及的并且計算效率很低．但是在模擬中全解析顆粒周圍的流場并準確捕捉顆粒的運動規(guī)律對于研究顆粒與流場及顆粒間相互作用非常重要．ohnon和Teduyar］利用全解析數(shù)值模擬方法模擬100個顆粒在三維通

于邊界處理部分考慮到Wu和附加邊界應(yīng)力方法非常適合于處理運動邊界本文球形曲面邊界處理均使用該方法詳見文獻［1，16顆粒碰撞部分的處理采用Glowinskietl［7］處理方法該方法假設(shè)當2顆粒間的距離小于給定的閥值時顆粒間將產(chǎn)生排斥作用7．數(shù)值方法驗單顆粒在無限大流場區(qū)域沉降是一個經(jīng)典問題有許多學者對此問題進行了研究本文將就此問題對模型進行驗算其計算域為12.5d×12.5d ρ10gcm3ρ771gcm3，流ν0009cm2sg981cms2．圖1給出不同網(wǎng)格分辨率情況下顆粒沉降速度隨時間變化的曲線圖 道中的運動．Panetl［10］模擬研究1204顆粒的流化運動規(guī)律．Kajishima［11］采用DNS

［17-18］結(jié)果進行對比在網(wǎng)格分辨率為dx=d/8d/16d/32(dx 為計算網(wǎng)格間距時的誤差分究2048

別為1.620.580.48但是dx=d32互作用．Yin和Koch［12］

數(shù)為1～20相應(yīng)的微觀結(jié)構(gòu)．Caoetal［13］用格子Boltzmann方法詳細分析2顆粒在不同初始狀態(tài)下顆粒相互作從上面的研究進展中可以看到顆粒懸浮運動中顆粒與流場相互作用起到重要作用并且對多顆粒懸浮運動的研究還很有限．因此本文詳細研究多顆粒懸浮運動的規(guī)律并試圖了解顆粒與流場的相互作用對顆粒運動產(chǎn)生的影響．本文采用格子oltmann方法及BF曲邊界處理格式研究三維大流場(全周期邊界)中多顆粒的懸浮運動過解析顆粒周圍的流場并且捕捉顆粒的運動規(guī)律詳細分析顆粒群運動過中相關(guān)統(tǒng)計量的規(guī)律． Boltzmann方法簡格子oltmann方法是以編簡單天然并行易于處理復(fù)雜邊界為特色的介觀方法數(shù)十年來格子oltmann方法在復(fù)雜流動等領(lǐng)域有了11］．本文采用LGK模型進行模擬速度模型是315詳細介紹見文獻1141．對

計算效率則選擇dx=d8或d16作為后面計算1不同網(wǎng)格大小情況下單顆粒沉降速度ig1Singleparticlesettlingvelocityversustimewithdifferentmeshsizes在討論多顆粒沉降的特性之前首先對雙顆粒沉降進行驗證計算域為6d6d60d，d/6cm網(wǎng)格分辨率為dx=d/6在初始時刻，23d3d27d3d25算例中，ρ10gcm3，ρ=114gcm3，ν= 中 大學學 第33001cm2sg981cms22給出2顆粒的運動軌跡及不同時刻顆粒的位置圖，(Drafting)，經(jīng)過一段時間后后面顆粒追趕上前面顆粒并碰撞(Kissing)2顆粒一起向下運動并發(fā)生翻轉(zhuǎn)Tumbling)，2個顆粒分開這種現(xiàn)111］．2兩顆粒的運動軌跡a及兩顆粒在不同時刻的位置圖(b)－(e)ig2Trajectoriesoftwo3可以看出兩表面顆粒之間的距離隨時間變化曲線與Apteetl［19］的結(jié)果比較符合2顆粒分離時刻相對于Glowinskietl［7］的結(jié)果有所延遲但是可以看出序可以很好地捕捉到DKT現(xiàn)象并很好地模擬顆粒運動及顆粒間相互作3兩顆粒之間的距離隨時間變化曲線ig3Distancebetweenthetwoparticlesversus數(shù)值模擬結(jié)果及分正如前面提到的多顆粒沉降過中的顆粒運動及顆粒流場相互作用是很重要的研究方向，顆粒運動規(guī)律背后的機制研究仍然不夠．本部分模擬研究144個顆粒和1152個顆粒的沉降過

給出顆粒平均沉降速度流場平均速度顆粒群平均各向速度顆粒群平均動能和顆粒群脈動速度等相應(yīng)統(tǒng)計值借此來研究多顆粒的運動規(guī)本節(jié)的1152等大小顆粒群模擬計算域為8d×32d流體網(wǎng)格分辨率為dx=d8，計算域的3個方向均為周期性邊界因此當顆粒穿過邊界后將從對面的邊界再次回到區(qū)域內(nèi)流體密度、顆粒直徑重力加速度依次為ρ10gcm－3f16cmg981cms2，度的比為1.4流體運動黏度為ν=0.01cm2/．初始時刻顆粒群在計算域中均勻分布2個相鄰顆粒間的距離為2圖4給出不同球面日網(wǎng)格點數(shù)(NL=29082)對顆粒沉降速度結(jié)果的影響圖中所顯示速度為顆粒群相對于流體的速度日分辨率的影響主要體現(xiàn)在顆粒群最大沉降速度上面，NL=82時顆粒群在受到各自尾渦的影響下加速過更長且受到的阻力更小(數(shù)值計算上)，故而最大速度變大但經(jīng)過各顆粒碰撞后，大體上都達到3cm/s速度．由于我們地關(guān)注顆粒群體系穩(wěn)定后的狀態(tài)由于計算資源有限所以我們采取NL=82對顆粒群系統(tǒng)進行了較長時間的運算后續(xù)結(jié)果不加說明情況下均為NL=82時所得．同時圖4也清晰地顯示顆粒群經(jīng)歷了最初的線性加速顆粒相互影響的非線性加速顆粒相互作用造成的過和振蕩過與之前的研究規(guī)律一致．圖4球面日網(wǎng)格點數(shù)NL對結(jié)果的影ig eswiththedifferentLagrangian圖5給出顆粒群xyz方向平均速度與合速度可以發(fā)現(xiàn):顆粒群主要在z方向(即沉降方向)上有較大速度和波動，Uz與合速度Up大致相等，在某種意義上來說xy方向分速度對合速度Up第2 查露等:三維顆粒群沉降的格子Boltzmann模 5顆粒xyz方向平均速UxUyUz與合速g．5Averagevelocitiesinx，y，andz從圖6可以看出顆粒群系統(tǒng)開始由勢能向動能轉(zhuǎn)化系統(tǒng)動能上升之后由于顆粒碰撞以及顆粒與流體相互作用等使得顆粒系統(tǒng)平均動能減少之后系統(tǒng)平均動能穩(wěn)定在1×10－9J左右由此也可以說明之后較長時間段內(nèi)顆粒群系統(tǒng)處于一個相對穩(wěn)定的狀態(tài)．6顆粒群平均動能隨時間變化關(guān)系ig6Variationoftheaveragekineticenergyoftheparticlegroupwithtime本節(jié)還對顆粒沉降運動進行了統(tǒng)計量的分析顆粒的脈動速度是用來反映顆粒速度偏離顆粒群平均速度的度本文用顆粒的離散速度的標準差來定義顆粒的脈動速度:v'x=

碰撞劇烈混合劇烈度高;碰撞之后顆粒群 y方向值但仍高于xy方向值．7顆粒群脈動速度ig7Fluctuatingvelocityofparticle84.9s—布圖1為小于vi－2(i為xyz)—間18vi22到17—( 別代表的是vi－2，vi+2025cm為間隔的16個均分小區(qū)間．系列1為z方向系列2為x方向系列3為y方向可以發(fā)現(xiàn)3個方向速度都呈現(xiàn)大多數(shù)顆粒速度降落在中間區(qū)域，擁有較大速度或較小速度的顆粒都處于少數(shù)大致呈現(xiàn)類似于分布的規(guī)律此時也說明顆粒呈現(xiàn)出均勻混合的理想氣體速度特性從另一側(cè)面說．N1N1(vx，p－vxv'∑ ∑ ( –)p

=

8xyz方向速度分布圖∑N–p (—v2 g．∑N–p (—v2v'

，－－ 顆粒群中顆粒分布不均勻性本文其中vxvyvz 擇用統(tǒng)計學中的變異系數(shù)來描述．這里簡單介xy和z方向上顆粒群的平均速度p為顆粒的 一下變異系數(shù)這是用來描述離散度的統(tǒng)計N從圖7可以看出顆粒群xy方向脈動自始

量．變異系數(shù)定義為Cv ∑∑Nsi(–i–n)/–2自終基本保持一致，xy方向同性．最初顆粒之 式中，ni為第i個子區(qū)域內(nèi)的顆粒數(shù)目(以顆－ 心點所在子區(qū)域為準)，N為子區(qū)域的數(shù)目， 用z值逐漸大于xy方向值表明此時顆粒之 子區(qū)域內(nèi)的平均顆粒數(shù)目．變異系數(shù)可以用來 中 大學學 第33示顆粒分布的非均勻性其值越大則表示非均勻性越強值為0時表示在子區(qū)域劃分的分辨率下圖90xy方向變異系數(shù)隨時間的變化規(guī)律可以發(fā)現(xiàn)xy方向走勢基本一致這也說明xy方向同性．在1s之后由于顆粒的相互作用諸如DKT過的影響顆粒開始向xy方向運動填補初始情況下邊界周邊的空地從而減小了非均勻性之后由于顆粒群的互相碰撞等強烈相互作用使得顆粒系統(tǒng)逐漸混合均勻趨向于理想原子運動狀態(tài)故而非均勻性進一步減小x、y變異系數(shù)值進一步減少．9x方向變異系數(shù)ig9Coefficientofvariationinx圖11顯示z方向的變異系數(shù)隨時間的變化關(guān)系初始階段的上下突變?yōu)檩^多顆粒同時穿越分割線造成并無多大實際意義由大體趨勢可以知道剛開始顆粒y向速度相差不大故相對位置變化不大v值變化不大之后由于顆粒的相互作用碰撞導致速度差異從而導致位置差異即

10y方向變異系數(shù)g．10Coefficientofvariationiny11z方向變異系數(shù)g．11CoefficientofvariationinzCv上升0.2左右三維多顆粒沉降顆粒流場相互作用圖12顯示不同時刻顆粒群的運動位置速度和中間截面流體速度黑色箭頭代表速度矢量圓球顏色代表各自的絕對速度標度尺顯示在圖13(a)的右下部分．t=0s時顆粒群和流體均處于靜止狀態(tài)．之后顆粒開始運動逐漸影響中間流12不同時刻顆粒及中間截面流體速度i12Fluidandparticlevelocityinthemiddlesectionatdifferent第2 查露等:三維顆粒群沉降的格子Boltzmann模 下弱流場上呈現(xiàn)的是紊亂的不規(guī)則的尾渦，顆粒之間發(fā)生相互作用．圖13(a)4.9s時顆粒群運動狀態(tài)俯視圖(b)4.9s時顆粒群運動矢量俯視ig13(a)Topviewoftheparticlegroupat4.9s;(b)motionvectorsat4.9s場截面t=0.7s時中間截面只受部分周圍顆粒對它的影響形成規(guī)律的一條條豎紋．之后顆粒群繼續(xù)發(fā)生相互作用越來越多的顆粒進入中間截面進而影響中間截面的流動狀況，t14s時有部分顆粒已近乎與中間流場速度保持一致．t=2.1s時顆粒群內(nèi)已不再像1.4s時還存在大量不等速度顆粒中間截面也是類似情況．顆粒群與流體之間的速度差異也越來越小從圖上來看可以大致認為t=2.1s時系統(tǒng)已處于速度空間的穩(wěn)定狀態(tài)．之后顆粒之間由于相互碰撞進一步在位置空間達到穩(wěn)定使整個顆粒群系統(tǒng)中顆粒具有類似于理想原子的運動規(guī)律表現(xiàn)為均勻而又混沌的狀態(tài)如圖13所示．圖14(a)和圖14(b)分別給出1152顆粒和114顆粒穩(wěn)定后某一時刻的渦量圖其中144顆粒算例初始情況除計算區(qū)域和各方向上顆粒數(shù)目縮短為1152顆粒的一半外其余均無變化．從圖中可看出顆粒周圍都出現(xiàn)較大的渦量值與an］結(jié)果吻合．還可看出在同樣初始體積分數(shù)下1000數(shù)量級顆粒群更容易出現(xiàn)大量顆粒群的團聚行為雖然100數(shù)量級的顆粒群也捕捉到了顆粒的KT現(xiàn)象和顆粒群的行為但規(guī)模不大多為34個顆粒．并且隨著顆粒數(shù)目的增加顆粒之間碰撞更易發(fā)生對周圍流場的擾動更頻繁顆粒對遠處的影響遠比少顆粒數(shù)目情況

圖14顆粒穩(wěn)態(tài)后中間截面渦量圖ig14Vorticitymapinmiddlesectionof1152particles(a)and114particles(b)結(jié)本文采用LB-EBF方法對三維懸浮顆粒群進行全解析直接數(shù)值模擬最多模擬1152顆粒(1000數(shù)量級)，對三維顆粒群的沉降過及其與顆粒群xy方向脈動自始自終基本保持一致，xy方向同性．最初較短時間內(nèi)各向同性，之后顆粒相互作用z值大于xy方向值．顆粒群系統(tǒng)處于穩(wěn)定狀態(tài)時速度分布大致呈現(xiàn)類似于分布的規(guī)律中間多兩頭少，顆粒呈現(xiàn)出均勻混合的理想氣體速度特性．xy方向變異系數(shù)走勢基本一致，xy方向同性，CvxCvy經(jīng)歷先減少后增大再減小的過．Cvz大致經(jīng)歷先增大再振蕩的過．1000數(shù)量級比100數(shù)量級顆粒群更容易出現(xiàn)大量顆粒群的團聚行為顆粒對遠處的影響遠比少顆粒數(shù)目情況下弱流場上呈現(xiàn)的是紊亂［1］AidunCK，ClausenJRLatticeboltzmannmethodforcomplexlowsAnnualReviewofFluidMechanics，2010，42:4392［2］AvciB，WriggersPADEMFEMcouplingapproachforthedirectnumericalsimulationof3Dparticulateflows［JJournalofAppliedMechanics，2012，79(1)10901107［3］FengZG，MichaelidesEETheimmersedboundary 中 大學學 第33Boltzmannmethodforsolvingfluidparticles betweenparticleclustersandparticleinducedturbulencproblem．JournalofComputationalPhysics，2004，195(2):6028．InternationalJournalofHeatandFluidFlow，2004，25(5):7218．［4FortesAF，JosephDD，LundgrenT．NonlinearYinXL，KochDLHinderedsettlingvelocityoffluidizationofbedsofsphericalparticle．JournalmicrostructureinsuspensionsofsolidsphereswithFluidMechanics，1987，177:467number．2．［5LomholtS，StenumB，MaxeyMRExperimentalCaoCS，ChenS，LiJ，etaSimulatingtheinctionsofthecouplingmethodforparticulateflows［JtwolysettlingsphericalparticlesinNewtonianfluidInternationalJournalofMultiphaseFlow，2002，28(2):latticeBoltzmannmethod［J．AppliedMathematicsComputation，2015，250:533［6FengJ，HuHH，JosephDDDirectsimulationofinitialvalueproblemsforthemotionofsolidbodiesinaudPart1Sedimentatio．JournalofFluid等流體動力學的格子Boltzmann方法M．:科學技術(shù)，2002:16KandhaiD，KoponenA，HoekstraA，eta1994，261:95aspectsof3DlatticeBGK:boundaries，accuracy，anda［7GlowinskiR，PanTW，HeslaTI，etlAfictitiousfastrelaxationmetho．JournalofComputationalapproachtothedirectnumericalsimulation 1999，150(2):4821viscousflowpastmovingrigidbodies:applicationWuJ，AidunCKSimulating3Ddeformableparticulateflo．JournalofComputationalsuspensionsusinglatticeBoltzmannmethodwith169(2):363externalboundary［J．InternationalJournal［8GaoH，LiH，WangLLatticeBoltzmannsimulationNumericalMethodsinFluids，2010，62(7):765