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

董輝，任旭云

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

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

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

1 DEM-CFD耦合法數學模型

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

1.1 液相控制方程

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

(1)

動量守恒方程為：

(2)

能量守恒方程為：

(3)

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

1.2 顆粒相控制方程

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

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

(4)

(5)

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

1.3 CFD-DEM耦合

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

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

2 數值仿真模型及驗證

2.1 幾何模型及參數

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

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

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

表1 顆粒參數Table 1 Parameters of the particles

2.2 網格獨立性驗證

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

圖3 不同網格數量砂顆粒沉降高度曲線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。采用上述實驗工況模擬單顆粒鋼球沉降，并將模擬仿真的結果與實驗結果作對比，如圖5所示。圖5為鋼球降落的垂直高度變化曲線，從圖5中可以看出DEM-CFD耦合模擬仿真與實驗結果基本一致，說明采用的耦合方法具有可行性。

圖4 鋼球實驗示意圖Fig.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為單個砂顆粒沉降過程中高度和速度曲線，黑色方塊曲線為砂顆粒的沉降高度變化曲線，綠色圓形曲線為砂顆粒的沉降速度變化曲線。從圖中可以看出，砂顆粒釋放后，在重力的作用下加快，直至碰到容器底部發生反彈。由于砂顆粒在水中運動受到阻力的影響，最終會在容器底部靜止。砂顆粒在到達容器底部之前加速降落，速度方向與重力方向一致。當砂顆粒與容器底部發生碰撞時，砂顆粒的速度方向變成向上，達到反彈最高點后，砂顆粒的沉降速度變負值，即開始再次沉降，多次碰撞反彈后砂顆粒在容器底部靜止。圖7為單個砂顆粒沉降過程中顆粒受力和角速度曲線，黑色三角形曲線為砂顆粒沉降過程中受力變化曲線，紅色星形曲線為砂顆粒沉降過程中角速度變化曲線。從圖中可以看出，砂顆粒受力沿重力方向逐漸減小，當與容器底部發生碰撞時，砂顆粒的受力突然增大后逐漸減小，再次發生碰撞時受力又會增加直至砂顆粒在容器底部靜止。砂顆粒在發生碰撞前角速度為0，在0.18 s左右發生碰撞，瞬時角速度增加到0.08 rad/s，說明顆粒與容器底部發生接觸后發生了旋轉，直到0.22 s砂顆粒在容器底部靜止時角速度變為0。

圖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，說明顆粒受力達到平衡，結合速度曲線也可看出水合物顆粒在0.2 s左右受力平衡進行勻速運動。水合物顆粒的角速度始終為0，說明顆粒在上升的過程中沒有發生旋轉。

圖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左右幾乎全部運動到容器底部，如圖10(c)所示。而水合物顆粒群全部運動到容器上部大約需要3.0 s，比砂顆粒群運動時間久，如圖10(f)所示。

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

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

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

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

4 結論

本文采用DEM-CFD耦合法，對砂顆粒和水合物顆粒進行三維數值模擬，主要結論如下：

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

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

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

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