基于群體平衡理論的管內(nèi)水合物沉積特性數(shù)值模擬

宋光春，李玉星，王武昌，姚淑鵬，魏 丁，閆 斌

（1.山東省油氣儲運安全省級重點實驗室 中國石油大學（華東），山東 青島 266580；

2.泰能天然氣有限公司，山東 青島 266580；3.煙臺新奧燃氣發(fā)展有限公司，山東 煙臺 264002）

水合物顆粒在油氣輸送管線內(nèi)的著床沉積是導致管線堵塞的重要原因［1］，因此，研究管道內(nèi)水合物的沉積特性對深水流動安全保障具有十分重要的意義［2-3］。目前，環(huán)道實驗是研究管道內(nèi)水合物沉積特性的主要方法［4］。近年來，隨著計算流體力學（CFD）的發(fā)展，數(shù)值模擬已成為替代環(huán)道實驗進行水合物沉積特性研究的有效手段。

Jassim等［5］使用CFD的方法研究了水合物在輸氣管線內(nèi)的沉積特性并建立了相應的沉積位置計算模型，但該模型并未實現(xiàn)水合物顆粒與流場的直接耦合。魏丁等［6］使用Fluent軟件耦合固液兩相流并以此模擬水合物顆粒在管線內(nèi)的沉積，模擬結(jié)果表明，流速對水合物沉積過程影響顯著，當流速達到一定值時，水合物體積分數(shù)對水合物沉積過程的影響不再重要。采用同樣的方法，陳鵬等［7-10］研究了流速、水合物體積分數(shù)、水合物顆粒粒徑、水合物顆粒黏度、壁面條件及管徑等因素對水合物沉積層高度的影響。但上述所有模擬均未考慮水合物顆粒粒徑的變化。Fatnes等［11］通過引入水合物顆粒聚集粒徑模型［12］計算水合物顆粒粒徑的變化并使用ANSYS CFX模擬了彎管處水合物的沉積情況，但該模型取水合物顆粒間的聚集力為常數(shù)且并未考慮顆粒間的聚并效率。Balakin等［13］使用群體平衡模型與CFD相結(jié)合的方法模擬了水合物在管道內(nèi)的聚集和沉積特性，該方法較全面地考慮了水合物顆粒在流動過程中的聚集和破碎，模擬結(jié)果與工程實際最接近。因此，群體平衡模型-CFD方法是對管道內(nèi)水合物顆粒沉積特性進行數(shù)值模擬的有效方法。

本工作引入了基于水合物顆粒聚集動力學的群體平衡模型，通過碰撞頻率和聚并效率分析水合物顆粒間的聚集，通過破碎頻率和破碎后子顆粒的粒徑分布函數(shù)分析水合物顆粒的破碎。利用Fluent 14.5軟件，對固液兩相流模型和群體平衡模型進行耦合求解，進而模擬流速、水合物顆粒粒徑及水合物體積分數(shù)對管內(nèi)水合物顆粒沉積特性的影響。

1 數(shù)值模型

1.1 幾何模型

為了精確模擬水合物在水平圓管內(nèi)的沉積特性，參考文獻［14］所用實驗環(huán)道建立了三維水平圓管模型（見圖1），管道長度為3 m，內(nèi)徑為45.2 mm。對幾何模型進行六面體網(wǎng)格劃分，進口壁面以1 mm為步長進行劃分；在近壁面處，為了處理邊界層效應，采用8層網(wǎng)格加密，其余網(wǎng)格均以1 mm為步長進行劃分。本工作幾何模型總計劃分267 696個六面體網(wǎng)格，網(wǎng)格質(zhì)量0.913。

圖1 水平圓管幾何模型及網(wǎng)格劃分示例Fig.1 Three-dimensional geometric model and mesh structure of the horizontal pipe.

對網(wǎng)格數(shù)量進行獨立性檢驗，結(jié)果見圖2。從圖2可看出，當網(wǎng)格數(shù)分別為527 616和267 696個時，管流（流速2 m/s）充分發(fā)展段近壁面處速度梯度的相對偏差較?。?.3%），說明網(wǎng)格數(shù)量滿足獨立性要求。

圖2 網(wǎng)格獨立性檢驗Fig.2 Grid dependency test.

1.2 多相流模型及湍流模型

本工作在建模過程中采用了以下假設：流動過程等溫且忽略相間質(zhì)量傳遞，即不考慮水合物的生成及分解；水合物漿由水相和水合物顆粒相構(gòu)成，不考慮氣相和油相；水合物顆粒粒徑連續(xù)分布；不考慮水合物顆粒在管壁上的黏附過程。在上述假設下，本工作采用的物理模型主要包括多相流模型、湍流模型和群體平衡模型。其中，多相流模型采用歐拉-歐拉雙流體模型，由連續(xù)性方程（見式（1））、動量方程（見式（2））兩個控制方程和若干用于封閉方程組的本構(gòu)方程構(gòu)成。

式中，i為水相或水合物顆粒相；t為時間，s；ρ為密度，kg/m3；α為體積分數(shù)，%；d為拉普拉斯算子；u為速度矢量，m/s；p為壓力，Pa；τi為應力張量，Pa；Mi則為相間動量交換項，kg/（m·s）2。

水合物漿的液固耦合是多相流模型建立時考慮的重點，在模擬過程中，液固耦合通過相間動量交換實現(xiàn)。在計算Mi時，主要考慮相間拖曳力（Mdi）和湍流擴散力（Mti），并以此建立相間作用力模型以封閉多相流模型，見式（3）。

式中，ur為相間相對流速，m/s；μtm為湍流黏度，kg/（m·s）；σd為Planck擴散系數(shù)；αs為水合物體積分數(shù)，%；kls為動量傳遞系數(shù)，本工作采用Gidaspow模型［15］計算：當αs≤20%時，由Wen-Yu公式計算，見式（4）；當αs＞ 20%時，由Ergun公式計算，見式（5）。

式中，L為水合物顆粒粒徑，m；μl為水相動力黏度，kg/（m·s）；CD為曳力系數(shù)。

除了相間作用力模型，多相流模型的封閉還需確定水合物顆粒相的黏度μs，計算式見式（6）［16］：

式中，αl為水相體積分數(shù)，%；μm為混合相黏度，可用Roscoe-Brinkman方程［17-18］計算得到（見式（7））。

根據(jù)式（6）～（7）即可編制用戶定義函數(shù)（UDF），并以此計算水合物顆粒相的黏度μs。

湍流模型則采用標準k-ε模型，近壁面采用標準壁面函數(shù)進行處理。

1.3 群體平衡模型

根據(jù)本工作假設，不考慮水合物顆粒的生成分解及水合物顆粒的管壁黏附且認為水合物顆粒粒徑連續(xù)分布，因此可采用式（8）所示的群體平衡方程［19］：

式中，n（L，t）表示粒徑為L的水合物顆粒在t時刻的數(shù)量密度，m-3；β（L-L′，L′）表示粒徑分別為L-L′和L′的水合物顆粒的碰撞頻率，m3/s；a為兩水合物顆粒發(fā)生碰撞后的聚并效率；S（L′）表示粒徑為L′的水合物顆粒的破碎頻率，s-1；b（L|L′）則表示粒徑為L′的水合物顆粒破碎后產(chǎn)生粒徑為L水合物顆粒的概率。

基于水合物顆粒聚集動力學，對式（8）中關(guān)鍵參數(shù)的計算公式進行選取。

碰撞頻率方面，主要考慮由差速沉降和流動剪切造成的碰撞并取兩碰撞頻率之和作為水合物顆粒的實際碰撞頻率。其中，差速沉降碰撞頻率（βDS）采用式（9）［20］計算：

式中，V為沉降速率，m/s，可由式（10）計算得到。

流動剪切碰撞頻率方面，當水合物顆粒小于Kolmogorov尺度時，處于湍流耗散區(qū)，在這一區(qū)域內(nèi)水合物的聚集行為主要受渦內(nèi)局部剪切力的影響。此時，水合物顆粒的碰撞頻率可用式（11）［21］進行計算；當水合物顆粒大于Kolmogorov尺度時，處于湍流慣性區(qū)，在這一區(qū)域內(nèi)水合物被主流場牽引運動。此時，水合物顆粒的碰撞頻率可用式（12）［22］進行計算：

式中，G為流場局部的絕對速度梯度，s-1；u為水合物顆粒的平均速度，m/s。

采用曲線模型計算水合物顆粒間的聚并效率，由于連續(xù)相為水，水合物顆粒間不存在液橋力，故計算聚并效率時主要考慮范德華力與流動剪切力之比。聚并效率的計算式見式（13）［23］：

式中，H為表征范德華力大小的哈梅克常數(shù)，J；R為發(fā)生碰撞兩水合物顆粒的調(diào)和半徑，m。

計算破碎頻率時則主要考慮由流動剪切導致的破碎，計算式見式（14）：

式中，E和m為經(jīng)驗常數(shù)。

破碎后子顆粒粒徑的分布函數(shù)方面，本工作采用二元分布作為水合物顆粒破碎后的粒徑分布函數(shù)［24-25］，表達式見式（15）。

根據(jù)式（9）～（15）編制UDF并以此計算群體平衡模型中的關(guān)鍵參數(shù)。

1.4 模型求解

利用Fluent 14.5軟件，對固液兩相流模型和群體平衡模型進行聯(lián)合求解。管道入口邊界條件設定為速度入口，管道出口邊界條件設定為壓力出口，以二階迎風差分格式離散動量方程，壓力-速度耦合使用SIMPLE算法，群體平衡方程則使用離散法進行求解，瞬態(tài)模擬時間步長設置為0.1 s，殘差設置為10-5。

模型求解時所涉及的部分參數(shù)見表1。

表1 模型參數(shù)Table 1 Model parameter

2 結(jié)果與討論

2.1 模擬工況

本工作主要研究流速、水合物顆粒粒徑及水合物體積分數(shù)對管內(nèi)水合物顆粒沉積特性的影響，具體模擬工況見表2。

表2 模擬工況Table 2 Simulation condition

2.2 實驗驗證

流動壓降是表征管道內(nèi)固液兩相流流動行為和流動特性的重要參數(shù)。根據(jù)相關(guān)實驗［14］，本工作選取流動壓降數(shù)據(jù)對數(shù)值模擬結(jié)果的可行性和有效性進行實驗驗證。表3為單位管長壓降實驗值和模型模擬值的對比。由表3可知，兩者的變化趨勢相同且兩者間的相對誤差均小于18%。因此，本工作采用的數(shù)值模型和模擬方法能較好地模擬水合物顆粒在管道內(nèi)的流動行為和流動特性。

表3 實驗和模擬條件下單位管長壓降對比Table 3 Comparison of pressure gradient between experimental and simulation results

2.3 流速對沉積特性的影響

水合物顆粒的著床沉積主要取決于顆粒所受的慣性力、浮力和重力，對于模擬工況下的水合物，它的密度大于水，故當重力大于浮力和曳力時，水合物顆粒就會沉積，形成床層。由表1可知，本工作中水合物顆粒的最大填充率為0.55，即當管道內(nèi)水合物體積分數(shù)達到55%時，水合物就會沉積形成床層。工況1～5的模擬結(jié)果見表4。從表4可看出，當管內(nèi)水合物初始粒徑為5 mm、水合物體積分數(shù)為10%時，只有流速為0.05 m/s的工況1發(fā)生了水合物沉積，管內(nèi)局部水合物最高體積分數(shù)達到了55%。圖3為工況1管道縱截面水合物體積分數(shù)云圖。由圖3可知，水合物沉積主要發(fā)生在模擬管段的中后部，在距離模擬管段入口約2.7 m的地方沉積現(xiàn)象最為明顯。

表4 工況1～5的模擬結(jié)果Table 4 Simulation results of case 1-5

圖3 工況1管道縱截面水合物體積分數(shù)云圖Fig.3 Contour of hydrate volume fraction distribution in pipeline longitudinal section of case 1.

從表4還可發(fā)現(xiàn)，流速的增加對管內(nèi)水合物顆粒粒徑大小的影響不具有明顯的規(guī)律性，這可能是因為：一方面，由水合物顆粒聚集動力學可知，流動剪切是導致水合物顆粒發(fā)生碰撞并產(chǎn)生聚集的主要原因，因此流速越高，水合物顆粒就越容易發(fā)生碰撞聚集從而導致粒徑增大；另一方面，流動剪切同樣是導致水合物顆粒發(fā)生剪切破碎的主要原因，因此流速越高，水合物顆粒也越容易發(fā)生剪切破碎從而導致粒徑減小。流速對管道橫截面水合物沉積情況的影響見圖4。從圖4可看出，流速還會對管內(nèi)水合物分布產(chǎn)生影響，流速越高，水合物分布就越均勻，水合物顆粒間發(fā)生碰撞聚集從而導致粒徑增大的概率也就越低。綜上所述，模擬工況的低流速范圍內(nèi)，流速的增加對管內(nèi)水合物顆粒粒徑大小的影響并無明顯的規(guī)律性。

由圖4還可看出，在工況1～5距管道入口2.7 m處的管道橫截面上，水合物均呈不均勻分布的特征，其中管道下部水合物濃度較高，管道中部水合物濃度分布相對均勻，管道上部水合物濃度則相對較低。固液流型方面，工況1出現(xiàn)了水合物沉積，為固定床層流，工況2～5則為不均勻懸浮流。隨著流速的增大，管道下部水合物濃度逐漸降低，管內(nèi)固液流型從固定床層流（有沉積）逐漸轉(zhuǎn)變?yōu)椴痪鶆驊腋×鳎o沉積），整個管道橫截面上水合物的濃度梯度也逐漸減小，水合物分布漸趨均勻。這是因為在管內(nèi)水合物體積分數(shù)和顆粒粒徑大小相近的情況下，流速越高，流體對水合物顆粒的攜帶和分散能力就越強，因此水合物顆粒發(fā)生著床沉積的可能性也就越小。

圖4 流速對管道橫截面水合物沉積情況的影響Fig.4 Influence of flow rate on hydrate deposition in pipeline cross section.

工況1距管道入口2.7 m處管道橫截面水合物顆粒粒徑分布云圖見圖5。

圖5 工況1管道橫截面水合物顆粒粒徑分布云圖Fig.5 Counter of hydrate particle diameter distribution in pipeline cross section of case 1.

由圖5可知，發(fā)生沉積時，管道橫截面處的水合物顆粒粒徑分布與水合物分布保持了較好的一致性。管道下部水合物顆粒粒徑較大，管道上部水合物顆粒粒徑較小，管道中部水合物顆粒粒徑分布則較均勻。這除了因為大粒徑水合物顆粒更易發(fā)生著床沉積外，還與管道下部水合物濃度高，水合物顆粒更易發(fā)生碰撞聚集從而導致粒徑增大有關(guān)。

圖6為工況1模擬管道各部分的水合物顆粒粒徑分布直方圖。由圖6a可知，模擬開始時，管道內(nèi)的水合物顆粒大小均在35 μm以下。之后，隨著水合物顆粒在流動過程中不斷聚集和破碎，管道內(nèi)水合物顆粒的粒徑逐漸變化到5～350 μm之間，近似呈對數(shù)正態(tài)分布（見圖6b）。最后，由于流速較低，顆粒粒徑大于70 μm的水合物顆粒幾乎全部發(fā)生著床沉積，只有粒徑小于70 μm的水合物顆粒在管道出口流出（見圖6c）。

圖6 工況1水合物顆粒粒徑分布直方圖Fig.6 Hydrate particle size distribution histogram of case 1.

2.4 水合物顆粒粒徑對沉積特性的影響

工況1及工況6～9的模擬結(jié)果見表5，水合物顆粒初始粒徑對水合物起始沉積位置的影響見圖7。從表5可看出，當管內(nèi)水合物體積分數(shù)為10%、管道入口流速為0.05 m/s時，水合物初始粒徑在5～200 μm之間的5組工況全部發(fā)生了水合物著床沉積。從圖7可看出，不同于工況1的是，工況6～9中水合物顆粒的起始沉積位置均距離管道入口較近，且水合物顆粒的初始粒徑越大，水合物顆粒的起始沉積位置距管道入口越近。此外，從表5還可看出，隨著模擬時水合物顆粒初始粒徑的增大，模擬得到的水合物顆粒平均粒徑和最大粒徑也逐漸增大。

表5 工況1及工況6～9的模擬結(jié)果表Table 5 Simulation results of case 1 and case 6-9

圖7 水合物顆粒初始粒徑對水合物起始沉積位置的影響Fig.7 Influence of initial hydrate particle diameter on initial hydrate deposition position.

圖8為工況1及工況6～9距管道入口2.7 m處管道橫截面水合物體積分數(shù)云圖。由圖8可知，由于出現(xiàn)水合物顆粒的著床沉積，上述5組工況對應的固液流型均為固定床層流。其中，管道下部水合物濃度較高形成固定床，中部水合物濃度分布相對均勻，上部水合物濃度則相對較低，幾乎無水合物顆粒存在（體積分數(shù)均低于0.004%）。由表5可知，隨著水合物顆粒初始粒徑的增大，模擬過程中管內(nèi)水合物顆粒的平均粒徑及最大粒徑也逐漸增大，因此管內(nèi)水合物的沉積情況漸趨嚴重，具體表現(xiàn)為管道下部高濃度水合物區(qū)和管道上部低濃度水合物區(qū)的范圍逐漸增大，管道橫截面上水合物的濃度梯度也越來越大。因此，在其他條件相同時，管內(nèi)水合物顆粒粒徑越大，管道橫截面處水合物分布越不均勻，管內(nèi)水合物沉積現(xiàn)象也越明顯。

值得注意的是，由圖8可看出，從工況1到工況6～9，管道橫截面處水合物沉積層的高度并未出現(xiàn)明顯變化。這是因為，工況1和工況6～9的床層高度為距管道入口2.7 m的管道橫截面處的床層高度。而由圖7及前文分析可知，不同工況下水合物在管道內(nèi)的起始沉積位置不同，同一工況下不同管道長度處水合物沉積層高度也不相同。因此，距管道入口2.7 m的管道橫截面并不是所有工況下水合物在管道內(nèi)沉積層高度最大的位置，故圖8中各個工況的床層高度并不能反映各個工況下管道內(nèi)水合物沉積的嚴重程度。不同工況下管道內(nèi)水合物沉積的嚴重程度應當從管道內(nèi)水合物沉積層的長度和沿程高度（即水合物沉積總量）來判斷，而不能單獨從某一截面上水合物沉積層的高度進行準確判斷。綜上所述，即使從工況1到工況6～9，管內(nèi)水合物的沉積情況越來越嚴重，但在管道某一長度的橫截面處（比如距入口2.7 m處），各工況的水合物沉積層高度仍可能相差不大。

圖8 顆粒粒徑對管道橫截面水合物沉積情況的影響Fig.8 Influence of hydrate particle diameter on hydrate deposition in pipeline cross section.

2.5 水合物體積分數(shù)對沉積特性的影響

工況1及工況10～13的模擬結(jié)果見表6。從表6可看出，當管內(nèi)水合物顆粒初始粒徑為5 μm、管道入口流速為0.05 m/s時，水合物體積分數(shù)在10%～50%之間的5組工況全部發(fā)生了水合物著床沉積。同工況6～9類似，工況10～13中水合物顆粒的起始沉積位置均距離管道入口較近，且水合物體積分數(shù)越大，水合物顆粒的起始沉積位置距管道入口越近。從表6還可看出，隨著模擬時水合物體積分數(shù)的增大，模擬得到的水合物顆粒的平均粒徑和最大粒徑也逐漸增大。

表6 工況1及工況10～13的模擬結(jié)果Table 6 Simulation results of case 1 and case 10-13

圖9為工況1及工況10～13距管道入口2.7 m處管道橫截面水合物體積分數(shù)云圖。由圖9可知，由于出現(xiàn)水合物顆粒的著床沉積，上述5組工況對應的固液流型均為固定床層流。其中，管道下部水合物濃度較高并形成固定床，管道中部水合物濃度分布相對均勻，管道上部水合物濃度則相對較低。

圖9 水合物體積分數(shù)對管道橫截面水合物沉積情況的影響Fig.9 Influence of hydrate volume fraction on hydrate deposition in pipeline cross section.

表7為工況1及工況10～13的管內(nèi)水合物分布。如表7和圖9所示，隨著管內(nèi)水合物體積分數(shù)的增大，管內(nèi)高水合物濃度（50%～55%）區(qū)所占的比例逐漸增大。因此，在其他條件相同時，管內(nèi)水合物體積分數(shù)越大，水合物沉積現(xiàn)象越明顯。

表7 工況1及工況10～13的管內(nèi)水合物分布Table 7 Hydrate volume fraction distribution of case 1 and case 10-13

3 結(jié)論

1）引入基于水合物顆粒聚集動力學的群體平衡模型，并將該模型與固液兩相流模型相結(jié)合，通過Fluent軟件進行聯(lián)合求解可實現(xiàn)對管道內(nèi)水合物顆粒沉積特性的精確模擬。

2）著床沉積發(fā)生時，管道橫截面處的水合物分布和粒徑分布有較好的一致性。其中，管道下部水合物濃度高、粒徑大，管道上部水合物濃度低、粒徑小，管道中部水合物的濃度和粒徑分布則較為均勻。

3）流速的增大會減小管道橫截面上水合物的濃度梯度，使水合物分布漸趨均勻，從而減弱水合物的著床沉積。

4）管內(nèi)水合物顆粒的初始粒徑越大、體積分數(shù)越高，水合物顆粒平均粒徑和最大粒徑的模擬值越大，水合物顆粒的起始沉積位置距管道入口越近，管內(nèi)水合物的沉積情況越嚴重。

符 號 說 明

a 聚并效率

b 子顆粒粒徑分布函數(shù)

CD曳力系數(shù)

E 經(jīng)驗常數(shù)

G 絕對速度梯度，s-1

g 重力加速度，m/s2

H 哈梅克常數(shù)，J

kls動量傳遞系數(shù)，

L，L′ 水合物顆粒粒徑，m

M 相間動量交換，kg/（m·s）2

m 經(jīng)驗常數(shù)

n 數(shù)量密度，m-3

p 壓力，Pa

R 調(diào)和半徑，m

S 破碎頻率，s-1

t 時間，s

u 速度矢量，m/s

u 平均速度，m/s

ur相間相對流速，m/s

V 沉降速率，m/s

α 體積分數(shù)

β 碰撞頻率，m3/s

ρ 密度，kg/m3

σdPlanck 擴散系數(shù)

τ 應力張量

μ 動力黏度，Pa·s

下標

1，2 分別表示不同的水合物顆粒

DS 差速沉降

di 相間拖曳力

i 水相或水合物相

l 水相

m 混合相

ST 層流剪切

s 水合物相

TH 湍流耗散區(qū)

TG 湍流慣性區(qū)

t 湍流

ti 湍流擴散力

［1］ Sloan E D. Natural gas hydrates in fl ow assurance［M］.New York：Elsevier Science Ltd，2010：1-36.

［2］ Song Guangchun，Li Yuxing，Wang Wuchang，et al. Experimental study of hydrate dissociation in oil-dominated systems using a high-pressure visual cell［J］.J Nat Gas Sci Eng，2017，45：26-37.

［3］ 宋光春，李玉星，王武昌，等. 油氣輸送管線水合物沉積研究進展［J］.化工進展，2017，36（9）：3164-3176.

［4］ Song Guangchun，Li Yuxing，Wang Wuchang，et al. Investigation of hydrate plugging in natural gas+diesel oil+water systems using a high-pressure fl ow loop［J］.Chem Eng Sci，2017，158：480-489.

［5］ Jassim E，Abdi M A，Muzychka Y. A new approach to investigate hydrate deposition in gas-dominated fl owlines［J］.J Nat Gas Sci Eng，2010，2（4）：163-177.

［6］ 魏丁，王武昌，李玉星，等. 管道 CCl3F 水合物漿流動特性的數(shù)值模擬［J］.油氣儲運，2016，35（8）：828-832.

［7］ 陳鵬，劉福旺，李玉星，等. 水合物漿液流動特性數(shù)值模擬［J］.油氣儲運，2014，33（2）：160-164.

［8］ 劉海紅. 水合物顆粒受力及聚結(jié)特性研究［D］.青島：中國石油大學（華東），2014.

［9］ 李瑩玉. 天然氣水合物顆粒在流場中的聚集行為研究［D］.天津：中國民航大學，2013.

［10］ 王曉婭. 管輸水合物漿液流動特性研究［D］.天津：中國民航大學，2015.

［11］ Fatnes E D. Numerical simulations of the flow and plugging behavior of hydrate particles［D］.Bergen：University of Bergen，2010.

［12］ Camargo R，Palermo T. Rheological properties of hydrate suspensions in an asphaltenic crude oil［C］//International Conference on Gas Hydrates，Yokohama：ICGH4，2002.

［13］ Balakin B V，Lo S，Kosinski P，et al. Modelling agglomeration and deposition of gas hydrates in industrial pipelines with combined CFD-PBM technique［J］.Chem Eng Sci，2016，153：45-57.

［14］ Balakin B V，Pedersen H，Kilinc Z，et al. Turbulent fl ow of freon R11 hydrate slurry［J］.J Pet Sci Eng，2010，70（3/4）：177-182.

［15］ Ding Jianmin，Gidaspow D. A bubbling fluidization model using kinetic theory of granular fl ow［J］.AIChE J，1990，36（4）：523-538.

［16］ 王繼紅. 冰漿的管道輸送熱流動特性［D］.大連：大連理工大學，2013.

［17］ Pabst W. Fundamental considerations on suspension rheology［J］.P R Soc A，2004，48（1）：6-13.

［18］ 趙鵬飛，王武昌，李玉星，等. 管道內(nèi)水合物漿流動的數(shù)值模型［J］.油氣儲運，2016，35（3）：272-277.

［19］ Hulburt H M，Katz S. Some problems in particle technology：A statistical mechanical formulation［J］.Chem Eng Sci，1964，19（8）：555-574.

［20］ Camp T R，Stein P C. Velocity gradients and internal work in fl uid motion［J］.J Bsn Soc Civ Eng，1943，30（4）：219-237.

［21］ Saffman P G，Turner J S. On the collision of drops in turbulent clouds［J］.J Fluid Mech，1956，1：16-30.

［22］ Abrahamson J. Collision rates of small particles in a vigorously turbulent fluid［J］.Chem Eng Sci，1975，30（11）：1371-1379.

［23］ van de Ven T G M，Mason S G. The microrheology of colloidal dispersions Ⅶ. Orthokinetic doublet formation of spheres［J］.Colloid Polym Sci，1977，255（5）：468-479.

［24］ Zhang Jianjun，Li Xiaoyan. Modeling particle-size distribution dynamics in a fl occulation system［J］.AIChE J，2003，49（7）：1870-1882.

［25］ 李振亮. 基于群體平衡的活性污泥絮凝動力學［D］.重慶：重慶大學，2014.

［26］ Li Xiaoyan，Logan B E. Collision frequencies between fractal aggregates and small particles in a turbulently sheared fluid［J］.Environ Sci Technol，1997，31（4）：1237-1242.