董輝，任旭云

(四川信息職業(yè)技術(shù)學(xué)院 智能控制學(xué)院，四川 廣元 628017)

海底非成巖水合物流化開(kāi)采過(guò)程中，形成了由水合物顆粒、砂顆粒和海水三相組成的水合物漿體[1-2],水合物漿體通過(guò)密閉管道如流化床反應(yīng)器、螺旋輸送機(jī)等[3-4]從海底輸送至海上工作平臺(tái)。為減少輸送的能源消耗和提高工作的可靠性，需要對(duì)采掘破碎后的水合物漿體中的泥砂進(jìn)行海底預(yù)分離。顆粒運(yùn)動(dòng)對(duì)海底非成巖水合物流化開(kāi)采具有十分重要的影響,因此研究水合物漿體海底泥砂預(yù)分離時(shí)，水合物和泥砂的顆粒運(yùn)動(dòng)是一個(gè)關(guān)鍵問(wèn)題。許多學(xué)者對(duì)顆粒運(yùn)動(dòng)現(xiàn)象進(jìn)行了研究，目前研究方法多是采用實(shí)驗(yàn)法和數(shù)值模擬法[5-7]，通過(guò)高速攝像技術(shù)研究顆粒在流體中的運(yùn)動(dòng)軌跡、沉降速度等。彭德其等[8]通過(guò)搭建豎管內(nèi)單顆粒自由沉降實(shí)驗(yàn)平臺(tái)，對(duì)顆粒沉降過(guò)程進(jìn)行了研究，得到了顆粒運(yùn)動(dòng)規(guī)律。Zhu等[9]開(kāi)發(fā)了一種無(wú)線球形探測(cè)器，通過(guò)該儀器測(cè)得球形顆粒在三維顆粒流中的平移和旋轉(zhuǎn)運(yùn)動(dòng)。粒子圖像測(cè)速技術(shù)(particle image velocimetry，PIV)[10]已應(yīng)用于顆粒的研究中，有許多學(xué)者通過(guò)該技術(shù)從圖像中提取并分析了流場(chǎng)、顆粒位置。通過(guò)實(shí)驗(yàn)的方法可以直觀地得到顆粒運(yùn)動(dòng)軌跡及位置數(shù)據(jù)，但是實(shí)驗(yàn)的成本較高，并且實(shí)驗(yàn)過(guò)程比較繁瑣，相比而言數(shù)值模擬方法具有成本低廉、分析條件靈活等優(yōu)點(diǎn)?；诖?，許多學(xué)者通過(guò)數(shù)值模擬的方法對(duì)顆粒運(yùn)動(dòng)進(jìn)行研究。Ma等[11]通過(guò)離散元素法(discrete element method，DEM)和計(jì)算流體力學(xué)(computational fluid dynamics，CFD)的耦合方法(DEM-CFD耦合法)和實(shí)驗(yàn)對(duì)比研究了單顆粒在旋流場(chǎng)中運(yùn)動(dòng)及碰撞特性，表明流場(chǎng)的旋轉(zhuǎn)速度和單個(gè)粒徑在顆粒運(yùn)動(dòng)行為中起著重要作用。Chen等[12]采用DEM-CFD耦合法研究了垂直管內(nèi)水力收集粗顆粒的顆粒運(yùn)動(dòng)及周圍流場(chǎng)。Hu等[13]利用格子玻爾茲曼方法研究了顆粒團(tuán)與單個(gè)顆粒在沉降過(guò)程中的相互作用、顆粒初始距離和顆粒數(shù)對(duì)顆粒團(tuán)運(yùn)動(dòng)和分布的影響。Cheng等[14]應(yīng)用DEM-CFD耦合法來(lái)模擬單通道泵中的固液流動(dòng)。

綜上所述，目前對(duì)顆粒運(yùn)動(dòng)規(guī)律的研究主要是針對(duì)顆粒沉降進(jìn)行速度與運(yùn)動(dòng)軌跡分析，但微觀上顆粒運(yùn)動(dòng)規(guī)律還研究的不夠深入，并且針對(duì)水合物顆粒這種密度小于水的顆粒運(yùn)動(dòng)研究的還比較欠缺。本文采用數(shù)值模擬耦合的方法，對(duì)砂顆粒和水合物顆粒分別進(jìn)行了單顆粒和顆粒群研究分析，表明了DEM-CFD耦合法能夠模擬大量顆粒運(yùn)動(dòng)的復(fù)雜流動(dòng)，并且能從微觀尺度分析顆粒受力和運(yùn)動(dòng)情況，為海底水合物藏開(kāi)采時(shí)水合物顆粒和砂顆粒運(yùn)動(dòng)的研究提供了參考，對(duì)不同密度的顆粒運(yùn)動(dòng)研究提供了一種可靠、有效的方法。

1 DEM-CFD耦合法數(shù)學(xué)模型

在DEM-CFD耦合仿真中，顆粒相用DEM方法進(jìn)行求解，液相用CFD方法求解，并通過(guò)固液兩相的相互作用實(shí)現(xiàn)耦合。

1.1 液相控制方程

水合物漿體多相流中的液相可視為黏性不可壓縮的流體，采用CFD求解，液相的質(zhì)量守恒方程為：

(1)

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

(2)

能量守恒方程為：

(3)

式中：ρf為液相密度，kg/m3；uf為液相流動(dòng)的速度，m/s；μf為液相的動(dòng)力黏度，N·s/m2；p為作用在液相微元體上的壓力，Pa；F為作用在流體微元體上的體積力，N；E為單位質(zhì)量液相的內(nèi)能，J；Ei為液相內(nèi)能，J；Ef為液相勢(shì)能，J。

1.2 顆粒相控制方程

水合物漿體多相流流動(dòng)過(guò)程中，砂和水合物顆粒均視為獨(dú)立的固相顆粒材料，可以采用DEM法進(jìn)行分析求解。DEM法主要應(yīng)用牛頓第二定律解析顆粒的運(yùn)動(dòng)，即可以描述每個(gè)顆粒的位移、速度和加速度。

本文中假定顆粒受到重力、顆粒-顆粒接觸力、顆粒-壁面接觸力、顆粒-流體相互作用力等。顆粒在t時(shí)刻的控制方程[15]為:

(4)

(5)

式中：m為顆粒的質(zhì)量，kg；u為顆粒的速度，m/s；Fp-p為顆粒間接觸力，N；Fp-w為顆粒與壁面產(chǎn)生的接觸力，N；Fp-f為顆粒與流體產(chǎn)生的接觸力，N；I為顆粒的轉(zhuǎn)動(dòng)慣量，kg·m2；ω為顆粒的角速度，rad/s；T為顆粒在質(zhì)心處受到的合外力矩，N·m。

1.3 CFD-DEM耦合

CFD與DEM耦合過(guò)程為：在每個(gè)時(shí)間步長(zhǎng)中，首先采用CFD求解計(jì)算域內(nèi)的流體，直至迭代計(jì)算收斂，然后根據(jù)顆粒所在的流體單元內(nèi)儲(chǔ)存的流體的物理參數(shù)(如速度、壓力等)，采用DEM計(jì)算出每個(gè)顆粒所受的力，并在相應(yīng)時(shí)間步長(zhǎng)內(nèi)進(jìn)行一次(或若干次)迭代計(jì)算。顆粒離散元仿真計(jì)算結(jié)束后，再次進(jìn)行流體CFD仿真計(jì)算。將所有顆粒的動(dòng)量匯作用在相應(yīng)的流體網(wǎng)格單元上，用于計(jì)算流體與顆粒之間相互能量的轉(zhuǎn)換，最終實(shí)現(xiàn)DEM-CFD雙向耦合，如圖1所示。

圖1 DEM-CFD耦合過(guò)程Fig.1 DEM-CFD coupling method

2 數(shù)值仿真模型及驗(yàn)證

2.1 幾何模型及參數(shù)

本文采用的是三維封閉容器，計(jì)算域?yàn)?00 mm×200 mm×200 mm的方形，如圖2(a)所示。重力方向?yàn)閦=-9.81 m/s2。模擬仿真流體域采用六面體結(jié)構(gòu)網(wǎng)格，如圖2(b)所示。

圖2 仿真計(jì)算域Fig.2 Simulation domain

基于DEM-CFD耦合仿真中，幾何模型為封閉的方形容器，內(nèi)部充滿靜止的水。壁面采用無(wú)滑移壁面邊界。顆粒與顆粒、顆粒與壁面的碰撞采用Hertz-Mindlin無(wú)滑動(dòng)接觸模型。仿真中顆粒參數(shù)如表1所示。

表1 顆粒參數(shù)Table 1 Parameters of the particles

2.2 網(wǎng)格獨(dú)立性驗(yàn)證

模擬計(jì)算采用六面體結(jié)構(gòu)網(wǎng)格，為驗(yàn)證仿真模型的網(wǎng)格獨(dú)立性，采用上述相同的物理參數(shù)和邊界條件對(duì)5種不同網(wǎng)格數(shù)量的模型進(jìn)行砂顆粒數(shù)值模擬，網(wǎng)格數(shù)量分別為80 256個(gè)、154 776個(gè)、204 511個(gè)、253 167個(gè)、310 823個(gè)，得到計(jì)算結(jié)果如圖3所示。從圖3中可以看出，網(wǎng)格數(shù)量對(duì)仿真結(jié)果有一定的影響，砂顆粒沉降曲線隨網(wǎng)格數(shù)量增加到一定程度后計(jì)算結(jié)果趨向穩(wěn)定。網(wǎng)格數(shù)量越多仿真計(jì)算所需要的時(shí)間越長(zhǎng)，因此綜合考慮仿真時(shí)間和仿真結(jié)果，選用網(wǎng)格數(shù)量為204 511個(gè)的模型可滿足網(wǎng)格獨(dú)立性。

圖3 不同網(wǎng)格數(shù)量砂顆粒沉降高度曲線Fig.3 Curves of particle settlement height with different grid numbers

2.3 仿真驗(yàn)證

為驗(yàn)證DEM-CFD耦合仿真的正確性，Li的實(shí)驗(yàn)工況是將直徑為9.5 mm的鋼球從靜止的甘油水溶液中釋放，如圖4所示[16]。鋼球的密度為7 780 kg/m3，泊松比為0.33，楊氏模量為200 GPa。采用上述實(shí)驗(yàn)工況模擬單顆粒鋼球沉降，并將模擬仿真的結(jié)果與實(shí)驗(yàn)結(jié)果作對(duì)比，如圖5所示。圖5為鋼球降落的垂直高度變化曲線，從圖5中可以看出DEM-CFD耦合模擬仿真與實(shí)驗(yàn)結(jié)果基本一致，說(shuō)明采用的耦合方法具有可行性。

圖4 鋼球?qū)嶒?yàn)示意圖Fig.4 Schematic of the steel sphere experiment

圖5 DEM-CFD 耦合仿真與實(shí)驗(yàn)比較Fig.5 Comparison of the experiment and DEM-CFD coupling

3 顆粒運(yùn)動(dòng)分析

3.1 砂顆粒沉降分析

對(duì)砂顆粒進(jìn)行耦合仿真，砂顆粒的粒徑為2 mm，將單個(gè)砂顆粒從高度為50 mm處釋放。圖6為單個(gè)砂顆粒沉降過(guò)程中高度和速度曲線，黑色方塊曲線為砂顆粒的沉降高度變化曲線，綠色圓形曲線為砂顆粒的沉降速度變化曲線。從圖中可以看出，砂顆粒釋放后，在重力的作用下加快，直至碰到容器底部發(fā)生反彈。由于砂顆粒在水中運(yùn)動(dòng)受到阻力的影響，最終會(huì)在容器底部靜止。砂顆粒在到達(dá)容器底部之前加速降落，速度方向與重力方向一致。當(dāng)砂顆粒與容器底部發(fā)生碰撞時(shí)，砂顆粒的速度方向變成向上，達(dá)到反彈最高點(diǎn)后，砂顆粒的沉降速度變負(fù)值，即開(kāi)始再次沉降，多次碰撞反彈后砂顆粒在容器底部靜止。圖7為單個(gè)砂顆粒沉降過(guò)程中顆粒受力和角速度曲線，黑色三角形曲線為砂顆粒沉降過(guò)程中受力變化曲線，紅色星形曲線為砂顆粒沉降過(guò)程中角速度變化曲線。從圖中可以看出，砂顆粒受力沿重力方向逐漸減小，當(dāng)與容器底部發(fā)生碰撞時(shí)，砂顆粒的受力突然增大后逐漸減小，再次發(fā)生碰撞時(shí)受力又會(huì)增加直至砂顆粒在容器底部靜止。砂顆粒在發(fā)生碰撞前角速度為0，在0.18 s左右發(fā)生碰撞，瞬時(shí)角速度增加到0.08 rad/s，說(shuō)明顆粒與容器底部發(fā)生接觸后發(fā)生了旋轉(zhuǎn)，直到0.22 s砂顆粒在容器底部靜止時(shí)角速度變?yōu)?。

圖6 砂顆粒沉降過(guò)程中高度和速度曲線Fig. 6 Curves of the particle height and velocity during sand particle settlement

圖7 砂顆粒沉降過(guò)程中受力和角速度曲線Fig. 7 Curves of the particle force and angular velocity during sand particle settlement

3.2水合物顆粒上升分析

對(duì)密度比水輕的水合物顆粒進(jìn)行耦合仿真，球形水合物顆粒的粒徑為2 mm，同樣將單個(gè)水合物顆粒從高度為50 mm處釋放。圖8為單個(gè)水合物顆粒上升過(guò)程高度和速度曲線，不考慮顆粒的融化，黑色方塊曲線為水合物顆粒的上升高度變化曲線，綠色圓形曲線為水合物顆粒的上升速度變化曲線。從圖中可以看出，水合物顆粒釋放后，由于水合物的密度比水小，顆粒在浮力的作用下上升。水合物顆粒從靜止釋放，在0.2 s左右有加速的過(guò)程，達(dá)到平衡后水合物顆粒勻速向上運(yùn)動(dòng)。圖9為單個(gè)水合物顆粒上升過(guò)程中顆粒受力和角速度變化曲線，黑色三角形曲線為水合物顆粒上升過(guò)程中受力變化曲線，紅色星形曲線為水合物顆粒上升過(guò)程中角速度變化曲線。從圖中可以看出，水合物顆粒沿重力反方向受力由5.8×10-6N逐漸減小，在0.2 s左右減小為0，說(shuō)明顆粒受力達(dá)到平衡，結(jié)合速度曲線也可看出水合物顆粒在0.2 s左右受力平衡進(jìn)行勻速運(yùn)動(dòng)。水合物顆粒的角速度始終為0，說(shuō)明顆粒在上升的過(guò)程中沒(méi)有發(fā)生旋轉(zhuǎn)。

圖8 水合物顆粒上升過(guò)程中高度和速度曲線Fig.8 Curves of the particle height and velocity during hydrate floating

圖9 水合物顆粒上升過(guò)程中受力和角速度曲線Fig.9 Curves of the particle force and angular velocity during hydrate floating

3.3 顆粒群分析

將密閉的容器中充滿混合均勻的砂顆粒和水合物顆粒，砂顆粒和水合物顆粒各有10 000個(gè)。顆粒群在重力、浮力和液體阻力共同作用下運(yùn)動(dòng)。圖10所示為不同時(shí)刻顆粒群的分布情況，圖10(a)為初始時(shí)刻砂顆粒和水合物顆粒均勻地分布在封閉容器中。當(dāng)顆粒群開(kāi)始釋放，由于砂顆粒的密度比水的密度大，會(huì)向容器底部運(yùn)動(dòng)，水合物顆粒密度比水的密度小，會(huì)向容器上部運(yùn)動(dòng)，如圖10(b)所示。砂顆粒群在0.6 s左右?guī)缀跞窟\(yùn)動(dòng)到容器底部，如圖10(c)所示。而水合物顆粒群全部運(yùn)動(dòng)到容器上部大約需要3.0 s，比砂顆粒群運(yùn)動(dòng)時(shí)間久，如圖10(f)所示。

圖10 不同時(shí)刻顆粒群分布情況Fig.10 Distribution of particles at different times

圖11為顆粒群的平均質(zhì)心位置變化曲線，圖12為顆粒群的平均速度變化曲線。初始時(shí)刻，砂顆粒群和水合物顆粒群均勻地分布在容器內(nèi)部，所以砂顆粒群和水合物顆粒群平均質(zhì)心都處在容器的中間位置。隨著顆粒群的釋放，砂顆粒群快速向容器底部沉降，水合物顆粒群緩慢向容器上部上升。砂顆粒群質(zhì)心位置的變化率比水合物顆粒的變化率大很多。從圖12也可看出相似的規(guī)律，砂顆粒群開(kāi)始運(yùn)動(dòng)時(shí)，平均速度突然增加，負(fù)值表示速度方向與重力方向一致。當(dāng)砂顆粒群逐漸運(yùn)動(dòng)到容器底部時(shí)，平均速度減小為0。相反，水合物顆粒群的平均速度增加的比較緩慢。導(dǎo)致砂顆粒群和水合物顆粒群變化率相差較大的原因是砂顆粒和水的密度差為1 600 kg/m3，而水合物顆粒和水的密度差僅有100 kg/m3左右。

圖11 顆粒群平均質(zhì)心位置曲線Fig.11 Curves of the average centroid location of particles

圖12 顆粒群平均速度曲線Fig.12 Curves of the average velocity of particles

4 結(jié)論

本文采用DEM-CFD耦合法，對(duì)砂顆粒和水合物顆粒進(jìn)行三維數(shù)值模擬，主要結(jié)論如下：

(1)單個(gè)砂顆粒在重力的作用下加速降落，直至碰到容器底部，當(dāng)與容器底部發(fā)生碰撞時(shí)砂顆粒的速度方向突然發(fā)生改變，受力突然增大，角速度瞬間增加，顆粒發(fā)生旋轉(zhuǎn)。

(2)由于水合物顆粒密度比水小，單顆粒在浮力的作用下上升，上升過(guò)程中在很短的時(shí)間加速，達(dá)到平衡后勻速上升。角速度始終為0，說(shuō)明水合物顆粒在上升的過(guò)程中沒(méi)有發(fā)生旋轉(zhuǎn)。

(3)砂顆粒和水合物顆粒群在重力、浮力和液體阻力共同作用下運(yùn)動(dòng)，砂顆粒群快速向容器底部沉降，水合物顆粒群緩慢向容器上部上升，其原因是砂顆粒與水的密度差大于水合物顆粒與水的密度差。

通過(guò)DEM-CFD耦合法可以模擬大量顆粒運(yùn)動(dòng)的復(fù)雜流動(dòng)，并且能從微觀尺度分析顆粒受力和運(yùn)動(dòng)情況，為水合物顆粒和砂顆粒運(yùn)動(dòng)的研究提供了一種可靠、有效的方法。本文在采用該方法進(jìn)行水合物和砂顆粒群共同運(yùn)動(dòng)的研究時(shí)，未分析顆粒組分分?jǐn)?shù)、顆粒粒徑以及液體密度對(duì)顆粒群運(yùn)動(dòng)的影響，這些不足之處有待繼續(xù)深入研究。