基于DEM-CFD耦合方法的水合物與砂顆粒運動分析

董輝，任旭云

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

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

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

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

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

1.1 液相控制方程

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

(1)

動量守恒方程為：

(2)

能量守恒方程為：

(3)

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

1.2 顆粒相控制方程

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

本文中假定顆粒受到重力、顆粒-顆粒接觸力、顆粒-壁面接觸力、顆粒-流體相互作用力等。顆粒在t時刻的控制方程[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)動慣量，kg·m2；ω為顆粒的角速度，rad/s；T為顆粒在質(zhì)心處受到的合外力矩，N·m。

1.3 CFD-DEM耦合

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

圖1 DEM-CFD耦合過程Fig.1 DEM-CFD coupling method

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

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

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

圖2 仿真計算域Fig.2 Simulation domain

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

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

2.2 網(wǎng)格獨立性驗證

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

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

2.3 仿真驗證

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

圖4 鋼球?qū)嶒炇疽鈭DFig.4 Schematic of the steel sphere experiment

圖5 DEM-CFD 耦合仿真與實驗比較Fig.5 Comparison of the experiment and DEM-CFD coupling

3 顆粒運動分析

3.1 砂顆粒沉降分析

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

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

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

3.2水合物顆粒上升分析

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

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

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

3.3 顆粒群分析

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

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

圖11為顆粒群的平均質(zhì)心位置變化曲線，圖12為顆粒群的平均速度變化曲線。初始時刻，砂顆粒群和水合物顆粒群均勻地分布在容器內(nèi)部，所以砂顆粒群和水合物顆粒群平均質(zhì)心都處在容器的中間位置。隨著顆粒群的釋放，砂顆粒群快速向容器底部沉降，水合物顆粒群緩慢向容器上部上升。砂顆粒群質(zhì)心位置的變化率比水合物顆粒的變化率大很多。從圖12也可看出相似的規(guī)律，砂顆粒群開始運動時，平均速度突然增加，負值表示速度方向與重力方向一致。當砂顆粒群逐漸運動到容器底部時，平均速度減小為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耦合法，對砂顆粒和水合物顆粒進行三維數(shù)值模擬，主要結(jié)論如下：

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

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

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

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