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

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

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

引言

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

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

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

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

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

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

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

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

取管道內(nèi)流體與顆粒群為研究對(duì)象,建立圓管內(nèi)流體與顆粒群耦合的兩相流模型,如圖1 所示.該模型共分2 個(gè)域,顆粒群域Ωp和流體域Ωf.流體域入口采用速度邊界Γin,出口采用壓力邊界Γout,速度和壓力用uΓin,f和pΓout,f表示.管道內(nèi)壁面邊界(流體域與管道域耦合界面)Γint假定為無(wú)滑移,用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]

流體動(dòng)量守恒方程為

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

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

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

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

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

1.2.1 顆粒與流體作用力分析

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

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

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

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

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

取任一顆粒i為研究對(duì)象,其運(yùn)動(dòng)由平移和旋轉(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 顆粒運(yùn)動(dòng)軌跡

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

顆粒在t+Δt時(shí)刻的運(yùn)動(dòng)位移為

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

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

1.3 顆粒碰撞模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中Q為碰撞前顆粒動(dòng)能.

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

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

1.4 耦合分析計(jì)算方法

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

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

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

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

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

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

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

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

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

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

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

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

圖6 顆粒碰撞搜索示意圖:(a)變網(wǎng)格建立示意圖;(b)顆粒運(yùn)動(dòng)軌跡穿過(guò)空間搜索網(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)顆粒碰撞列表建立

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

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

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

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

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

則顆粒j運(yùn)動(dòng)軌跡的參數(shù)方程為

式中φ 為參數(shù).

H平面的點(diǎn)法式方程為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.4 顆粒群碰撞求解流程

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

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

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

3 CFD-DEM 耦合分析及顆粒碰撞搜索算法驗(yàn)證

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

3.1 兩個(gè)顆粒碰撞算例

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

如圖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顆粒勻速運(yùn)動(dòng).工況2:A,B顆粒變速運(yùn)動(dòng),沿x方向加速度為-10 m/s2,沿y方向加速為-10 m/s2.根據(jù)顆粒碰撞理論和解析解,得到A,B顆粒的碰撞速度、碰撞位置和發(fā)生碰撞的時(shí)間、碰撞力如表1所示.

表1 兩顆粒運(yùn)動(dòng)及碰撞的理論和數(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 取計(jì)算時(shí)間0.1 s,計(jì)算時(shí)間步為10-4s,10-5s,10-6s,10-7s 對(duì)顆粒勻速、變速運(yùn)動(dòng)進(jìn)行數(shù)值模擬,得到碰撞力、碰撞速度、碰撞位置等結(jié)果一并列入表1.模擬所用計(jì)算機(jī)為英特爾Core i7-9700 @ 3.00 GHz 八核,主顯卡AMD Radeon RX 550 Series,內(nèi)存16 GB.

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

圖10 兩顆粒勻速運(yùn)動(dòng)誤差分析圖Fig.10 Error analysis diagram of two particles in uniform motion

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

圖11 兩顆粒變速運(yùn)動(dòng)誤差分析圖Fig.11 Error analysis diagram of variable speed movement of two particles

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

3.2 多個(gè)顆粒碰撞算例

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

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

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

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

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

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

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

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

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

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

圖15 多顆粒運(yùn)動(dòng)軌跡與碰撞示圖Fig.15 Movement trajectories and collisions of more than one particle

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

3.3 CFD-DEM 耦合分析算例

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

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

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

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

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

圖16 顆粒運(yùn)動(dòng)軌跡Fig.16 Multiple particle trajectories in fluid

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

4 結(jié)論

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

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

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