顆粒群碰撞搜索及CFD-DEM 耦合分域求解的推進算法研究1)

劉巨保 王 明,2) 王雪飛 姚利明 楊 明,3) 岳欠杯

*(東北石油大學(xué)機械科學(xué)與工程學(xué)院,黑龍江大慶 163318)

?(新加坡南洋理工大學(xué)機械與航天工程學(xué)院,新加坡 639798)

引言

在自然界和工業(yè)領(lǐng)域,流體中的顆粒碰撞普遍存在[1].如泥沙流中的砂粒碰撞與沉積[2],石油開采攜砂壓裂液中的顆粒碰撞與管壁沖蝕等.這些顆粒群與流體運移的兩相流動中,流體輸運的湍流效應(yīng)和顆粒運動的碰撞相互影響[3].在顆粒體積分?jǐn)?shù)低至4.0×10-4情況下[4],顆粒間的碰撞效應(yīng)也不應(yīng)被忽略.通常采用數(shù)值模擬手段來量化整個流場以及單個顆粒的動力學(xué)[5].

固液兩相流動中,顆粒間接觸和碰撞十分劇烈、復(fù)雜[6],成為影響兩相流動行為的關(guān)鍵因素[7],也是數(shù)值模擬研究的熱點和難點.固液兩相流的數(shù)值模型主要有歐拉-歐拉雙流體模型(TFM)[8-9]和歐拉-拉格朗日離散顆粒模型(CFD-DEM)[10]兩大類.前者將顆粒與流體視為互相滲透的擬流體或擬連續(xù)介質(zhì),在歐拉坐標(biāo)系研究;雙流體模型可以全面考慮顆粒相的輸運特性,適合進行大規(guī)模的工程問題計算,但忽略了顆粒離散特性,對顆粒流動行為的預(yù)測會與實際過程產(chǎn)生偏差[11],無法得到顆粒的運動軌跡[12].后者采用歐拉-拉格朗日坐標(biāo)系研究兩相流動,即流體作為連續(xù)介質(zhì),在歐拉坐標(biāo)系研究,顆粒作為離散相,在拉格朗日坐標(biāo)系基于牛頓運動定律求解顆粒動力學(xué)方程,跟蹤顆粒在流場中的復(fù)雜運動.歐拉-拉格朗日模型對流體與顆粒耦合計算的模擬精確程度,取決于能否準(zhǔn)確地刻畫顆粒間的碰撞,而顆粒碰撞搜索算法成為解決這一難題的關(guān)鍵.

前人對顆粒碰撞搜索算法進行了大量研究,如分割單元算法也稱網(wǎng)格單元算法[13-15],將計算域劃分成均勻的立方體單元,單元長度為顆粒直徑的整數(shù)倍.計算時間步結(jié)束時刻根據(jù)顆粒當(dāng)前位置被分配到某個單元,只有同一單元或相鄰單元內(nèi)顆?？赡馨l(fā)生碰撞.一個顆粒最多占據(jù)8 個單元,當(dāng)其成為目標(biāo)顆粒時,只對其占據(jù)的單元內(nèi)其他顆粒進行碰撞檢測.該算法遺漏了那些時間步結(jié)束前,離開或穿過立方體單元,可能與目標(biāo)顆粒發(fā)生碰撞的顆粒.Allen等[16]提出的鄰居列表算法,規(guī)定了一個閾值距離來確定顆粒間相互作用的影響半徑.碰撞搜索只對列表內(nèi)的顆粒對進行搜索,如果兩個顆粒之間的距離大于閾值距離,其就不可能在下一個時間步中發(fā)生碰撞.所以,鄰居列表需要幾個時間步更新一次.該算法構(gòu)建的搜索列表內(nèi)顆粒眾多,在高濃度顆粒存在的應(yīng)用中,每對顆粒進行一次判定計算代價過高.Vemuri 等[17]提出的八叉樹算法,建立一個能夠容納任意方向上最大顆粒的立方體,將其作為劃分域的基本單元.在對小顆粒進行搜索碰撞時,基本單元將會被多次劃分八叉樹,每個顆粒根據(jù)其大小與單元的關(guān)系,配置到對應(yīng)單元中.在計算域內(nèi)有多種尺寸顆粒時,需多次劃分基本單元,計算過程繁瑣.

Baraff[14]提出了邊界盒算法,使用軸向?qū)R的邊界盒來確定顆粒之間的碰撞.如果兩個顆粒相撞,其在x,y和z維度上的正交投影將會重疊.3 個列表包含對軸的投影坐標(biāo),每個維度對應(yīng)一個列表.在下一個時間步驟,需要再次對之前排序的列表進行排序.Li 等[18]提出的時間相干碰撞搜索算法(SMB),是將包圍盒AABB 投影在笛卡爾坐標(biāo)軸上,根據(jù)投影間距遞增排序建立初步碰撞搜索列表,通過對列表內(nèi)所有間距搜索檢查確定是否重疊,建立碰撞列表.該算法適用于大量移動AABB 的系統(tǒng),并支持對象的動態(tài)插入和刪除,不受物體大小的限制.

上面幾種方法對顆粒碰撞的判定,都采用較小的時間步長進行推進計算,根據(jù)計算時間步結(jié)束時顆粒位置,判斷顆粒是否與目標(biāo)顆粒碰撞.若時間步選取過大,會遺漏計算時間步內(nèi)的顆粒碰撞效應(yīng)[19],無法對顆粒碰撞導(dǎo)致的顆粒位置、速度、加速度等變化進行描述影響計算精度.而無限縮小顆粒計算時間步,又會造成推進次數(shù)的增加,影響計算效率.

此外,由于顆粒碰撞搜索方式的改變,顆粒與流體耦合方式也相應(yīng)發(fā)生變化.傳統(tǒng)CFD-DEM 耦合方法,流體計算時間步長是顆粒計算時間步長的整數(shù)倍(Δtf=n·Δtp)[20].流體完成一個時間步Δtf,顆粒完成n個時間步Δtp,顆粒與流體才進行一次相間力的交換,導(dǎo)致顆粒在流體內(nèi)部受力及運動與實際誤差較大.

針對CFD-DEM 耦合計算時,顆粒碰撞搜索時間步選取過小導(dǎo)致計算效率低、選取過大導(dǎo)致碰撞顆粒漏判的問題.本工作通過對顆粒碰撞搜索算法的研究,使顆粒計算時間步的選取盡可能大、且不受顆粒碰撞時間限制;依據(jù)顆粒與流體耦合條件,自適應(yīng)調(diào)整流體計算時間步長,實時更新相間力.以期為顆粒-流體耦合計算提供一種高計算精度、高計算效率的計算方法.

1 流體與顆粒群耦合的數(shù)值分析方法

取管道內(nèi)流體與顆粒群為研究對象,建立圓管內(nèi)流體與顆粒群耦合的兩相流模型,如圖1 所示.該模型共分2 個域,顆粒群域Ωp和流體域Ωf.流體域入口采用速度邊界Γin,出口采用壓力邊界Γout,速度和壓力用uΓin,f和pΓout,f表示.管道內(nèi)壁面邊界(流體域與管道域耦合界面)Γint假定為無滑移,用uΓI,f=0表示.在圓管內(nèi)兩相流體系中任取一控制體dVc,該控制體內(nèi)流體和顆粒占據(jù)的體積分別為dVf和dVp.

圖1 流體與顆粒群耦合模型示意圖Fig.1 Schematic diagram of fluid and particle group coupling model

1.1 流體控制方程

流體質(zhì)量守恒方程為[21]

流體動量守恒方程為

式中,ρf為流體密度,uf為流體速度矢量,pf為流體壓力,g為重力加速度,αf為流體的體積分?jǐn)?shù);為流體黏性應(yīng)力張量,Sf為流體與顆粒動量交換力源項.

兩相流與單相流的主要區(qū)別在于,流體相中流體體積分?jǐn)?shù)αf和源項Sf的變化.流體體積分?jǐn)?shù)αf是流體單元中流體所占有的體積份額[22],Sf是顆粒與流體間動量交換力源項[23],如圖2 所示.

圖2 兩相流體計算單元示意圖Fig.2 Schematic diagram of two-phase fluid calculation unit

式中,VE為流體單元體積,αp為流體單元內(nèi)顆粒占有的體積分?jǐn)?shù),為流體單元內(nèi)第i個顆粒體積,為流體單元內(nèi)第i個顆粒受到流體作用力,np為流體單元內(nèi)顆粒個數(shù);FE和ME分別為流體單元內(nèi)所有顆粒與流體作用力的合力和合力矩,H為顆粒到流體單元質(zhì)心的距離.

1.2 顆粒受力分析及動力學(xué)方程

1.2.1 顆粒與流體作用力分析

顆粒被裹挾在流體內(nèi)部,跟隨流體運動,不穩(wěn)定的流態(tài)會對顆粒產(chǎn)生多種作用力,其合力Ff可表示為

式中,Fd為顆粒在流體中受到的阻力(曳力),Fp為壓力梯度力,Fvm為虛擬質(zhì)量力[24],Fsaff為Saffman 剪切升力[25],Fml為Magnus 旋轉(zhuǎn)升力[26],Fb為倍瑟特力[27].

由于顆粒直徑較小,忽略Saffman 剪切升力、Magnus 旋轉(zhuǎn)升力、壓力梯度力、虛擬質(zhì)量力、倍瑟特力對顆粒的影響[28].阻力(曳力)在流體作用于顆粒上的力中起主要作用[29].當(dāng)顆粒的運動速度大于流體的流速時,流體作用在顆粒上的力表現(xiàn)為阻力,當(dāng)流體的流速大于顆粒速度時,作用在顆粒上的力表現(xiàn)為曳力.計算公式為[30]

式中,Cd為曳力系數(shù),Rep為顆粒雷諾數(shù),Dp為顆粒直徑,μf為流體表觀黏度.

1.2.2 顆粒動力學(xué)方程

取任一顆粒i為研究對象,其運動由平移和旋轉(zhuǎn)組成,根據(jù)牛頓第二定律得[31]

式中,mi為顆粒i質(zhì)量,Fcn,ij和Fct,i j分別為顆粒i與顆粒j的法向和切向碰撞力,Fg,i為顆粒i重力,Ii為顆粒i的慣性矩,ωi為顆粒i的角速度,Mc為顆粒碰撞引起的力矩

1.2.3 顆粒運動軌跡

假設(shè)顆粒在任意t時刻,合力為F(t),速度為u(t),位置為(x(t),y(t),z(t)).此時顆粒的加速度為a(t)=F(t)/mi.當(dāng)顆粒運動Δt后,其瞬時速度為[32]

顆粒在t+Δt時刻的運動位移為

式中,γ 為插值系數(shù),γ=0.5 為中心差分.

在空間直角坐標(biāo)系下,顆粒在t+Δt時刻顆粒的位置為(x(t+Δt),y(t+Δt),z(t+Δt)),其計算公式為

1.3 顆粒碰撞模型

顆粒碰撞過程的假設(shè):(1) 顆粒是球形的剛體;(2)顆粒間只存在二體碰撞;(3)顆粒密度相同、直徑可以不同,顆粒群初始狀態(tài)為均勻分布.為了建立顆粒間碰撞的計算方法,不妨取顆粒群中任兩顆粒記為顆粒i和顆粒j,設(shè)顆粒i直徑為質(zhì)量為mi、顆粒j直徑為質(zhì)量為mj,兩顆粒質(zhì)量比為q=mi/mj.

兩顆粒發(fā)生碰撞滿足動量定理,如圖3 所示[33].

圖3 顆粒碰撞示意圖Fig.3 Schematic diagram of particle collision

式中,ui和為顆粒i碰撞前、后平動速度;ωi和為顆粒i碰撞前、后轉(zhuǎn)動速度;k為顆粒相對速度的法向單位向量;J為顆粒碰撞時的沖量;Ii和Ij為顆粒i與顆粒j的轉(zhuǎn)動慣量.

當(dāng)兩個顆粒碰撞后有滑移時,碰撞后平動速度和轉(zhuǎn)動速度為

顆粒在碰撞時停止滑移,碰撞后平動速度和轉(zhuǎn)動速度為

式中,λ 為摩擦系數(shù),e為顆粒碰撞恢復(fù)系數(shù),τ為顆粒相對速度的切向單位向量,uij為顆粒碰撞前相對速度,ui j,τ為顆粒碰撞前相對速度的切向分量.

顆粒碰撞時產(chǎn)生滑移滿足下列條件

式中,Jτ為顆粒碰撞沖量的切向分量.

根據(jù)彈性力學(xué)碰撞理論[34],兩顆粒碰撞時間tc為

根據(jù)沖量定理,可得兩顆粒碰撞力計算公式為

由式(15) 和式(16) 可知,兩顆粒碰撞過程滿足動量守恒,但碰撞過程能量發(fā)生耗散,不考慮顆粒旋轉(zhuǎn)時的能量損失,則顆粒碰撞耗散能量為[35]

計算時間步內(nèi)顆粒受到流體作用力的耗散能量為

由式(28)和式(29),得到顆粒i能量耗散率為

式中Q為碰撞前顆粒動能.

顆粒碰撞計算流程如圖4 所示,碰撞顆粒i與顆粒j相對速度uij通過式(25) 判斷顆粒是否發(fā)生滑移,若發(fā)生滑移則通過式(17)～式(20)計算顆粒碰撞后速度,若顆粒不發(fā)生滑移則通過式(21)～式(24)計算顆粒碰撞后速度.采用插值算法,將顆粒碰撞前相對速度uij與碰撞后相對速度進行加權(quán)插值得到代入式(26)和式(27)得到顆粒碰撞時間tc和碰撞力Fc,i j,然后將其代入式(9)和式(10)得到顆粒碰撞后的速度,并對碰撞后速度的設(shè)定和求得值進行比較判斷,若不滿足需修正碰撞后速度,重復(fù)式(26)、式(27)和式(9)、式(10)計算;否則,求得顆粒碰撞后速度值.其中,h為求解次數(shù),為速度收斂容差.

圖4 顆粒碰撞后速度計算流程Fig.4 Calculation flowchart of velocity for Particle collision

1.4 耦合分析計算方法

顆粒位置變化直接導(dǎo)致流體體積分?jǐn)?shù)和動量交換力源項的變化,因此,根據(jù)流體單元內(nèi)流體體積分?jǐn)?shù)和動量交換力源項的變化,建立流體與顆粒群耦合的收斂條件

CFD-DEM 耦合分析時,其耦合計算流程如圖5所示,CFD 計算采用SIMPLE 算法求解流體域連續(xù)性方程和動量方程,根據(jù)顆粒和流體的相對速度,計算顆粒受到的流體作用力Ff,傳遞給DEM 求解器,并啟動DEM 計算.DEM 采用牛頓第二運動定律求解顆粒運動方程,根據(jù)顆粒在Δtp的運動軌跡,搜索和判斷顆粒碰撞,若顆粒產(chǎn)生碰撞,需修正顆粒計算時間步長.根據(jù)顆粒在計算域內(nèi)位置和受力的變化計算得到流體的體積分?jǐn)?shù)αf和動量交換力源項Sf,若不滿足式(31)和式(32),傳遞給CFD 求解器,自動更新流體計算時間步長,并推進CFD 計算.

CFD 求解如圖5(a) 所示,對流體計算域進行初始化,設(shè)置流體初始計算時間步長Δtf,判斷流體是否滿足收斂條件,若滿足,CFD 求解結(jié)束.若不滿足,縮短流體時間步Δtf=kn· Δtf(k≤1 為修正系數(shù),n=1,2,3,··· 為流體時間步縮短次數(shù)),重新進行CFD 計算,直至滿足流體收斂條件,CFD 求解結(jié)束.

CFD 求解結(jié)束后,計算Δtf時刻顆粒受到的流體力Ff傳遞給DEM,并啟動DEM 求解.

DEM 求解如圖5(b)所示,設(shè)置顆粒計算初始時間步長Δtp=Δtf,判斷顆粒在Δtp時刻計算結(jié)果是否滿足式(31)、式(32) 顆粒與流體耦合收斂條件,若不滿足,縮短顆粒時間步=kn· Δtp,重新進行DEM 求解,直至滿足耦合收斂條件,DEM 求解結(jié)束,此時耦合分析的時間步為若滿足,根據(jù)顆粒在Δtp的運動軌跡,搜索和判斷顆粒是否產(chǎn)生碰撞,若無碰撞,=Δtp,DEM 求解結(jié)束.若存在碰撞,計算得到最先發(fā)生碰撞的時間≤Δtp,并根據(jù)的顆粒位置,判斷是否滿足耦合收斂條件.若不滿足,縮短顆粒求解時間重新進行DEM 求解,直至滿足耦合收斂條件,DEM 求解結(jié)束,此時耦合分析的時間步為(其中i=1 時,時間步內(nèi)第1 次搜索到碰撞,若滿足,需判斷碰撞累計時長是否滿足若不滿足,DEM 求解結(jié)束.若滿足,則以剩余時間作為新的顆粒時間步長重新進行顆粒計算,并搜索顆粒碰撞,根據(jù)的顆粒位置,判斷是否滿足耦合收斂條件,若不滿足,縮短顆粒求解時間重新進行DEM 求解,直至滿足耦合收斂條件,DEM 求解結(jié)束,此時耦合分析的時間步為(其中i=2 時,時間步內(nèi)第2次搜索到碰撞,若滿足耦合條件,再次判斷碰撞累計時長是否滿足＜Δtp(i=2 時,若不滿足DEM 求解結(jié)束.若滿足,重新以剩余時間進行顆粒計算直至碰撞累計時長＜Δtp不成立,或者不滿足耦合收斂條件,此時(時間步內(nèi)進行n次碰撞)DEM 求解結(jié)束.

圖5 CFD-DEM 耦合計算流程圖Fig.5 Flow chart of CFD-DEM coupling calculation

DEM 求解結(jié)束后,計算流體體積分?jǐn)?shù)αf和動量交換力源項Sf,并根據(jù)耦合分析時間步對顆粒和流體計算時間步進行修正Δtp=判斷流體計算時間tf+Δtf＜T是否成立,若成立,將αf和Sf傳遞給CFD,重新開始CFD 求解;若不成立,則耦合計算結(jié)束.

其中,顆粒碰撞的搜索采用的是逆向迭代算法,增大計算時間步長不會影響顆粒碰撞搜索的計算精度.顆粒不發(fā)生碰撞時,計算時間步長的增大,使顆粒運動位移增大,導(dǎo)致流體計算單元內(nèi)流體體積分?jǐn)?shù)αf和動量交換力源項Sf變化過快,造成顆粒-流體耦合計算誤差.式(31)和式(32)限制誤差的大小,保證顆粒-流體耦合計算誤差不隨計算時間步長變化.

2 顆粒碰撞搜索的逆向迭代算法

流體內(nèi)部的大量顆粒群隨流體進行不規(guī)則運動,顆粒與顆粒、顆粒與壁面存在碰撞,快速、準(zhǔn)確判斷顆粒間發(fā)生碰撞,是提高顆粒與流體耦合計算效率的有效方法.

2.1 目標(biāo)顆粒的變網(wǎng)格構(gòu)造

在數(shù)值模擬顆粒間碰撞時,在一個時間步內(nèi)顆粒只能與周邊顆粒發(fā)生碰撞[36],將目標(biāo)顆粒與周圍可能發(fā)生碰撞的顆粒組成一個搜索空間網(wǎng)格進行碰撞搜索.具體方法如下:將目標(biāo)顆粒在計算時間步的運動位移,作為空間網(wǎng)格對角線,建立搜索空間網(wǎng)格體,如圖6(a) 所示.由于計算域內(nèi)不同顆粒在Δt時間步內(nèi)運動的位移不同,所構(gòu)建的空間網(wǎng)格也會出現(xiàn)不同尺寸,即為變網(wǎng)格.

圖6 顆粒碰撞搜索示意圖:(a)變網(wǎng)格建立示意圖;(b)顆粒運動軌跡穿過空間搜索網(wǎng)格示意圖;(c)顆粒間距篩選示意圖Fig.6 Search diagram of particle collision:(a)Schematic diagram of variable grid establishment;(b)schematic diagram of particle trajectories passing through the spatial search grid;(c)schematic diagram of particle spacing screening

2.2 目標(biāo)顆粒碰撞列表建立

為了縮短顆粒碰撞搜索和判定時間,對目標(biāo)顆粒周圍顆粒篩選,保留與目標(biāo)顆?？赡馨l(fā)生碰撞的顆粒,建立碰撞搜索列表,具體步驟如下.

(1)在計算時間步內(nèi),篩選運動軌跡與空間網(wǎng)格體存在交集的顆粒,如圖6(a)所示.

設(shè)目標(biāo)顆粒在t時刻和t+Δt時刻的坐標(biāo)為(x0,y0,z0)和(x1,y1,z1),以此兩點連線為對角線,構(gòu)建目標(biāo)顆粒的碰撞搜索空間網(wǎng)格體.

對計算域內(nèi)任一顆粒j,從t到t+Δt時刻的運動軌跡與目標(biāo)顆粒所建空間網(wǎng)格體若有交集,則可能產(chǎn)生碰撞;否則無碰撞.顆粒運動軌跡穿過空間網(wǎng)格體必定與空間網(wǎng)格體對角面H或W相交,因此計算顆粒運動軌跡與對角面的交點,并判斷交點是否在空間網(wǎng)格體的范圍內(nèi)即可,如圖6(b)所示.

不妨取顆粒j運動軌跡經(jīng)過坐標(biāo)(xj,yj,zj),且方向向量為m=(Xj,Yj,Zj),在H平面經(jīng)過點(x0,y0,z0),法向量為n=(n1,n2,n3),運動軌跡與H平面交點o(x,y,z).

則顆粒j運動軌跡的參數(shù)方程為

式中φ 為參數(shù).

H平面的點法式方程為

其中,顆粒j運動軌跡與H平面的交點o(x,y,z)一定滿足式(33)和式(34),將式(33)代入式(34)得

將φ 代入式(33) 得到交點o(x,y,z) 坐標(biāo),若交點坐標(biāo)在空間網(wǎng)格體范圍內(nèi),則說明運動軌跡穿過空間網(wǎng)格體.若交點坐標(biāo)不在空間網(wǎng)格體范圍內(nèi),或式(35) 分母為0 (運動軌跡與H平面平行),則運動軌跡與H平面不相交.需用相同方法判斷運動軌跡與W平面是否相交,若運動軌跡與W平面相交,說明運動軌跡穿過空間網(wǎng)格體.

(2)依據(jù)顆粒與目標(biāo)顆粒運動軌跡間距離小于或等于兩顆粒半徑之和的條件,再次篩選可能碰撞顆粒,如圖6(c)所示.

目標(biāo)顆粒i運動軌跡方向向量為e=(Xi,Yi,Zi),顆粒j運動軌跡方向向量為b=(Xj,Yj,Zj).將e×b=c⊥得到向量e,b的公垂向量c⊥=(X,Y,Z).

取目標(biāo)顆粒i運動軌跡上任意一點E和顆粒j運動軌跡上任意一點B,得到向量eb,將向量eb在公垂向量c⊥方向上做投影,得到目標(biāo)顆粒i與顆粒j運動軌跡間最短距離d如下

在計算時間步內(nèi),若兩顆粒運動軌跡最近距離d≤則目標(biāo)顆粒i可能與顆粒j發(fā)生碰撞,顆粒j保留.若d＞則目標(biāo)顆粒i不會與顆粒j發(fā)生碰撞.

經(jīng)(1)(2)兩步篩選后,保留下來的顆粒形成目標(biāo)顆粒碰撞搜索列表.

2.3 碰撞顆粒的逆向搜索算法

在提高顆粒碰撞計算效率時,采用了較大的計算時間步Δtp后,為了解決顆粒碰撞漏判狀況,沿顆粒走過的運行軌跡采用逆向搜索方式,搜索判斷兩顆粒是否存在碰撞.若不存在,說明該時間步內(nèi)計算結(jié)果正確,可進行下一個時間步計算,若存在碰撞,確定發(fā)生碰撞時間,修正顆粒計算時間步.

如圖7 所示,假設(shè)在計算時間步Δtp內(nèi)存在t+Δtc時刻,顆粒i與顆粒j產(chǎn)生碰撞.根據(jù)顆粒運動軌跡計算方法,得到顆粒i與顆粒j在時間增量為Δtc時顆粒的位置坐標(biāo)為

圖7 顆粒碰撞時間求解示意圖Fig.7 Schematic diagram of particle collision time solution

顆粒i與顆粒j在Δtc處的距離和發(fā)生碰撞條件為

方程中顆粒i和顆粒j在t時刻的速度和坐標(biāo)均為已知,方程可化簡成一元四次方程

對方程求解得到Δtc1,Δtc2,Δtc3,Δtc4四個解,依據(jù)條件,判斷解的有效性.若存在真解,取最小值作為顆粒i與顆粒j碰撞時間步.若無真解,說明顆粒在計算時間步內(nèi)不發(fā)生碰撞.

在一個計算時間步內(nèi),假設(shè)目標(biāo)顆粒i與碰撞列表內(nèi)顆粒mi發(fā)生碰撞,其碰撞時間步為1,2,···,m).因只考慮顆粒二體碰撞,最早與目標(biāo)顆粒相遇的顆粒k才是發(fā)生碰撞的顆粒,則碰撞時間步取為

2.4 顆粒群碰撞求解流程

對于由n個顆粒組成的顆粒群,依據(jù)前述方法,可選擇n-1 個目標(biāo)顆粒,完成n-1 個變網(wǎng)格和碰撞搜索列表建立.假設(shè)目標(biāo)顆粒i對顆粒j進行搜索判斷后建立碰撞搜索列表,顆粒j作為目標(biāo)顆粒進行列表建立時,無需對顆粒i重新進行判斷,直接從總的碰撞列表內(nèi)映射即可.

將每個列表內(nèi)的最小碰撞時間步長,放入總表內(nèi),得到最小碰撞時間步,局部求解域的計算時間步逆向到最先發(fā)生碰撞時刻,執(zhí)行碰撞.碰撞計算完成后,碰撞顆粒重新更新碰撞列表,并計算最小碰撞時間步長,其他未發(fā)生碰撞顆粒,只需在原有碰撞時間步長基礎(chǔ)上減去已經(jīng)完成的碰撞時間步長,更新碰撞時間步總表.顆粒群碰撞求解流程如圖8 所示.

圖8 顆粒群碰撞搜索求解流程Fig.8 Particle swarm collision search solution process

3 CFD-DEM 耦合分析及顆粒碰撞搜索算法驗證

基于Fluent 和DEM 平臺,采用Delphi 編程二次開發(fā)了流體與顆粒群耦合動力學(xué)分析程序.

3.1 兩個顆粒碰撞算例

采用現(xiàn)有的DEM 軟件和本文提出的改進顆粒碰撞搜索算法MDEM 分別對兩個顆粒碰撞進行數(shù)值模擬.預(yù)期驗證本文方法的計算精度和計算效率.

如圖9 所示,取A,B兩顆粒的半徑R=5 mm,彈性模量為210 GPa,泊松比為0.27,A顆粒初始速度=2 m/s,=-2 m/s,B顆粒初始速度=2 m/s,A,B顆粒初始圓心間距100 mm.為了便于討論,記和分別為A,B顆粒碰撞后x,y方向速度,Fcn為A,B顆粒的碰撞力.工況1:A,B顆粒勻速運動.工況2:A,B顆粒變速運動,沿x方向加速度為-10 m/s2,沿y方向加速為-10 m/s2.根據(jù)顆粒碰撞理論和解析解,得到A,B顆粒的碰撞速度、碰撞位置和發(fā)生碰撞的時間、碰撞力如表1所示.

表1 兩顆粒運動及碰撞的理論和數(shù)值模擬結(jié)果Table 1 Theoretical and numerical simulation results of movement and collision of two particles

圖9 顆粒位置示意圖Fig.9 Schematic diagram of particle position

DEM 和MDEM 取計算時間0.1 s,計算時間步為10-4s,10-5s,10-6s,10-7s 對顆粒勻速、變速運動進行數(shù)值模擬,得到碰撞力、碰撞速度、碰撞位置等結(jié)果一并列入表1.模擬所用計算機為英特爾Core i7-9700 @ 3.00 GHz 八核,主顯卡AMD Radeon RX 550 Series,內(nèi)存16 GB.

由圖10 可知,在勻速工況下,采用DEM 方法,只有當(dāng)計算時間步在10-6s 及以下,得到的計算結(jié)果與理論解誤差才不大于2%,其中速度最大誤差為2%和1%,碰撞位置最大誤差均為0.2%,碰撞力誤差為1.61%和1.13%,發(fā)生碰撞時刻誤差為1.6%和1.5%,A2和A3計算耗時分別為24.73 s 和54.04 s.由此可見,DEM 方法的計算精度受計算時間步影響較大,計算效率隨著時間步的減小而顯著降低.采用MDEM 方法,不同計算時間步得到的計算結(jié)果與理論結(jié)果誤差均小于2%,其中速度最大誤差為1%,1%,0.5%,碰撞位置最大誤差均為0.22%,0.2%,0.2%,碰撞力誤差為0.85%,0.98%,0.88%,發(fā)生碰撞時刻誤差為0.9%,1.3%,0.09%,B1,B2,B3計算耗時分別為20.01 s,23.14 s,37.06 s.由此可見,本文提出的MEDM方法,計算精度不受計算時間步的影響,且隨著時間步的減小計算效率呈下降趨勢,其中B2計算時間步長比A2計算時間步長大一個量級,計算耗時縮短4.59 s (6.4%).B1計算時間步長比A2計算時間步長大兩個量級,計算耗時縮短4.72 s(19.1%).

圖10 兩顆粒勻速運動誤差分析圖Fig.10 Error analysis diagram of two particles in uniform motion

由圖11 可知,在變速工況下,采用DEM 方法,只有當(dāng)計算時間步在10-6s 及以下,得到的計算結(jié)果與理論解誤差均小于2%,其中速度最大誤差為1.930%和0.400%,碰撞位置最大誤差均為0.17%和0.38%,碰撞力誤差為1.62%和1.13%,發(fā)生碰撞時刻誤差為1.6%和1.4%,C2和C3計算耗時27.73 和58.21s.采用MDEM 方法,不同計算時間步得到的計算結(jié)果與理論結(jié)果誤差均小于2%,其中速度最大誤差為0.26%,0.26%,0.13%,碰撞位置最大誤差均為0.15%,0.45%,0.12%,碰撞力誤差為1.64%,1.32%,0.9%,發(fā)生碰撞時刻誤差為0.9%,0.9%,1.8%,D1,D2,D3計算耗時22.2 s,25.1 s,38.11 s.D2計算時間步長比C2計算時間步長大一個量級,計算耗時縮短2.63 s(9.48%),D1計算時間步長比C2計算時間步長大兩個量級,計算耗時縮短5.53 s(19.9%).由此可見,兩種方法在模擬變速運動時的結(jié)論和規(guī)律同勻速運動.

圖11 兩顆粒變速運動誤差分析圖Fig.11 Error analysis diagram of variable speed movement of two particles

綜上所述,采用本文方法對顆粒勻速、變速運動進行數(shù)值模擬,選取不同時間步均能得到近似的計算結(jié)果,當(dāng)計算時間步比DEM 大兩個量級時,勻速運動的計算效率提高19.1%,變速運動的計算效率提高19.9%,使時間步選取過小、計算效率低的問題得到有效解決.

3.2 多個顆粒碰撞算例

顆粒幾何形狀、材質(zhì)屬性同3.1 算例.顆粒初始位置如圖12 所示,水平設(shè)置6 個顆粒,顆粒間距20 mm,編號h1,h2,h3,h4,h5,h6,其中h1,h3,h5顆粒水平速度均為0 m/s,垂直向下速度均為2 m/s,h2,h4,h6顆粒水平和垂直速度均為2 m/s.豎直設(shè)置6 個顆粒,顆粒間距20 mm,編號v1,v2,v3,v4,v5,v6,其中v1,v3,v5顆粒水平速度均為2 m/s,垂直速度均為0 m/s,v2,v4,v6顆粒水平和垂直速度均為2 m/s,顆粒做勻速運動.

圖12 多顆粒碰撞初始位置示意圖Fig.12 Schematic diagram of initial position of multi-particle collision

采用理論分析,可得6 對顆粒h1–h2,h3–h4,h5–h6,v1–v2,v3–v4,v5–v6首次發(fā)生碰撞時刻為0.005 0 s,碰撞位置分別為(25,10),(65,10),(105,10),(10,25),(10,65),(10,105),后繼顆粒碰撞較復(fù)雜,無法得到理論解.在DEM 和MDEM 數(shù)值模擬中,分別選取計算時間步為10-5s,10-6s,10-7s 和10-4s,10-5s,10-6s,計算時間取0.1 s.

由 DEM 數(shù)值模擬結(jié)果可知,計算時間步為10-5s 得到23 次顆粒對碰撞,而計算時間步為10-6s和10-7s,得到25 次顆粒對碰撞,部分結(jié)果見表2.由圖13 可知,在10-6s 和10-7s 計算時間步得到的首次發(fā)生碰撞時刻、碰撞位置與理論解的誤差均不大于2%.計算耗時分別為60.31 s 和90.12 s.由此可見,計算時間步取10-5s 時,DEM 方法漏判了2 次顆粒碰撞,隨著計算時間步的減小,多顆粒碰撞數(shù)值模擬的計算效率降低.

圖13 DEM 顆粒碰撞時間、位置誤差分析圖Fig.13 DEM particle collision time and position error analysis diagram

由MDEM 數(shù)值模擬結(jié)果可知,計算時間步取10-4s,10-5s,10-6s 時,均得到25 次顆粒對碰撞,部分結(jié)果如表2 和表3 所示.由圖14 可知,不同計算時間步得到的首次發(fā)生碰撞時刻、碰撞位置與理論解的誤差均不大于2%.計算耗時分別為51.31 s,56.22 s,65.21 s.由此可見,在模擬多顆粒碰撞時,MEDM 方法的計算精度不受計算時間步大小影響,可采用較大的時間步來提高計算效率.

圖14 MDEM 顆粒碰撞時間、位置誤差分析圖Fig.14 MDEM particle collision time and position error analysis diagram

表2 多顆粒運動及碰撞的數(shù)值模擬結(jié)果(DEM)Table 2 Numerical simulation results of multi-particle motion and collision(DEM)

表3 多顆粒運動及碰撞的數(shù)值模擬結(jié)果(MDEM)Table 3 Numerical simulation results of multi-particle motion and collision(MDEM)

當(dāng)DEM 和MDEM 計算時間步均取10-6s 時,MDEM 比DEM 計算時間增加4.9 s,這是由于每個計算時間步內(nèi)的逆向搜索耗時所致.當(dāng)MDEM 的計算時間步擴大一個量級(10-5s) 時,比DEM 計算時間縮短4.09 s,即計算耗時降低了6.8%;當(dāng)MDEM 的時間步取10-4s 時,比DEM 計算時間縮短了9 s,即計算耗時降低了14.9%.其中h4,v4,h2,v2均發(fā)生5 次碰撞,運動軌跡如圖15 所示,虛線為MDEM,實線為DEM.

圖15 多顆粒運動軌跡與碰撞示圖Fig.15 Movement trajectories and collisions of more than one particle

綜上所述,采用本文方法選取10-6s,10-5s,10-4s 時間步對多個顆粒碰撞進行數(shù)值模擬,均能得到精確的計算結(jié)果,計算時間步選取可比DEM 方法大兩個量級,其計算耗時降低了14.9%.

3.3 CFD-DEM 耦合分析算例

采用CFD-DEM 和CFD-MDEM 方法分別對3.2算例中的多個顆粒與流體耦合進行數(shù)值模擬.計算流體域為300 mm×300 mm×300 mm 的矩形域,流體密度ρ=1 g/cm3、黏度μf=0.1 g/(cm·s),流體初始狀態(tài)為靜止.計算域采用規(guī)則的矩形網(wǎng)格離散,除上部邊界設(shè)置為開邊界外,其余壁面采用無滑移邊界條件.顆粒幾何參數(shù)、材質(zhì)、初始位置及計算參數(shù)等與算例3.2 完全相同,CFD-DEM 流體計算時間步長設(shè)置為顆粒計算時間步長的10 倍.CFD-MDEM 顆粒計算初始時間步長與流體計算初始時間步相等,顆粒與流體真實計算時間步在耦合計算過程中,根據(jù)耦合收斂條件實時調(diào)整和修正,時間步長設(shè)置如表4所示.

表4 顆粒與流體耦合分析的時間步長設(shè)置Table 4 Time step setting of particle-fluid coupling analysis

CFD-DEM 計算結(jié)果表明:在A1,A2,A3計算時間步得到顆粒碰撞次數(shù)分別為18 次、22 次、22次,與算例3.2 相比碰撞次數(shù)減少了3 次,原因在于流體的阻力使顆粒做減速運動,計算耗時分別為18 min,27 min,42 min,只有在顆粒計算時間步小于等于10-6s 時,CFD-DEM 耦合計算才能得到精確值.

CFD-MDEM 計算結(jié)果表明:在B1,B2,B3計算時間步,得到的顆粒碰撞次數(shù)均為22 次,計算耗時分別為22.5 min,24 min,31 min,顆粒計算初始時間步的改變不影響CFD-MDEM 耦合計算精度.

A2流體計算時間步比B3流體計算初始時間步大一個量級,A2流體求解器啟動104次,顆粒求解器啟動105次,B3流體與顆粒求解器均啟動100 383 次.B3比A2多啟動383 次,是由于MDEM 算法根據(jù)耦合收斂條件,修正計算初始時間步,得到真實計算時間步比計算初始時間步小,B3比A2計算多耗時4 min.B2流體計算初始時間步與A2流體計算時間步相等,B2顆粒計算初始時間步比A2顆粒計算時間步大一個量級,B2流體與顆粒求解器啟均動11 479 次,與A2流體求解次數(shù)相比求解器多啟動1479 次,與B3相比求解器多啟動次數(shù)由383 次增加到1479 次,是由于顆粒計算初始時間步的擴大,顆粒求解不滿足收斂條件次數(shù)增多造成.B2比A2計算少耗時3 min,降低11.1%.B1流體計算初始時間步比A2流體計算時間步大一個量級,顆粒計算初始時間步比A2顆粒計算時間步大兩個量級,B1流體與顆粒求解器均啟動3254 次,比A2計算少耗時4.5 s,降低16.7%.A2和B1顆粒h2,h3,h4和v2,v4,v6的運動軌跡基本一致,如圖16 所示.

圖16 顆粒運動軌跡Fig.16 Multiple particle trajectories in fluid

綜上所述,采用CFD-MDEM 方法對多個顆粒與流體進行耦合數(shù)值模擬,流體與顆粒計算初始時間步取10-4s,10-5s,10-6s 均能得到精確的計算結(jié)果,計算耗時比CFD-DEM 方法降低了16.7%.

4 結(jié)論

本文針對CFD-DEM 耦合計算中,顆粒時間步選取過小、計算效率低,時間步選取過大、碰撞顆粒漏判的問題,通過算法研究和算例驗證,取得如下結(jié)論和認(rèn)識.

(1)提出了顆粒碰撞搜索的改進算法,以目標(biāo)顆粒運動位移構(gòu)建可變空間搜索網(wǎng)格,篩選可能碰撞顆粒建立搜索列表,減少了顆粒碰撞的搜索時間;采用逆向搜索方式判斷顆粒碰撞,避免了時間步選取過大對顆粒碰撞漏判,確保了顆粒碰撞的準(zhǔn)確描述.通過兩個顆粒和多個顆粒的運動和碰撞算例表明,本方法計算時間步選取的大小對顆粒碰撞計算精度影響不大;兩個顆粒在勻速和變速工況下的數(shù)值模擬計算耗時比DEM 方法降低了19.1%和19.9%,多個顆粒的數(shù)值模擬計算耗時比DEM 方法降低了14.9%.

(2)建立了顆粒與流體計算時間步自動調(diào)整的耦合算法,顆粒與流體計算初始時間步選取可以相同,根據(jù)顆粒推進時間實時調(diào)整流體計算時間,實現(xiàn)流體計算時間步隨顆粒計算時間自動調(diào)整,提高了顆粒與流體耦合的計算效率.多個顆粒與流體耦合算例的數(shù)值模擬表明,該算法的計算精度受顆粒和流體的初始計算時間步長影響較小,可通過較大的顆粒時間來提高計算效率,該算法在求解算例的耗時比CFD-DEM 方法降低了16.7%.