基于LBM 水力裂縫中支撐劑顆粒團簇升阻力系數研究

張 濤，吳春燕，孫 堃，盧 聰，李 聰

(1.西南石油大學石油與天然氣工程學院，成都 610000；2.空氣動力學國家重點實驗室，綿陽 621000)

石油工程中，油氣藏儲層水力壓裂時壓裂液攜支撐劑在縫內輸送是一個狹窄裂縫中的固液兩相流過程[1-3]，此過程中支撐劑顆粒會聚集成團形成不同構型的顆粒團簇[4]，從而影響顆粒阻力的大小[5]。同時，裂縫壁面的存在會影響顆粒團簇受力，從而對顆粒的運移鋪置產生影響[6]。因此，明確顆粒團簇在窄縫內運移過程中的受力規律，對水力壓裂過程中支撐劑輸送問題的研究有重要意義。

格子玻爾茲曼方法(LBM, Lattice Boltzmann Method)在介觀尺度將流體離散為微團進行模擬計算。相對于物理實驗和基于分子動力學的微觀尺度模擬以及基于N-S 方程宏觀尺度模擬，LBM 沒有連續性假設且計算量小，對復雜流固邊界附近的流場刻畫較為精細[7]，容易獲得顆粒在流場中的受力信息。多位學者采用該方法研究流固多相流動過程。BEETSTRA 等[5]基于LBM 推導不同幾何形狀顆粒群的曳力系數公式，其研究表明顆粒群中任意單顆粒的曳力系數低于孤立的單顆粒曳力系數。GONG 等[8]采用LBM 研究顆粒簇的阻力系數與顆粒雷諾數在沉降過程中的關系，隨著顆粒雷諾數的增加，團簇中顆粒間相互作用減少。HILL等[9]基于LBM 研究得到曳力隨孔隙率及雷諾數變化的關系式。YIN 等[10]基于LBM 對多分散相系下顆粒曳力系數進行了研究，并提出了曳力系數公式。SARKAR 等[11]模擬了液體中的顆粒曳力，擴大了固液多分散系相修正因子的應用范圍。SOMMERFELD 等[12]應用LBM 確定不規則形狀顆粒流動阻力系數的平均標準差隨雷諾數增加而減小，表明顆粒形狀的不規則性對較高的雷諾數有更強的影響。王春雨等[13]模擬了非均勻多孔顆粒的升力隨達西數變化規律。楊佼[14]采用D2Q9模型，研究了雙顆粒、四顆粒在不同顆粒間距、不同顆粒排布方式和不同雷諾數下的繞流模式與受力特性。

從上述文獻看出，國內外學者基于LBM 對單顆粒及顆粒群的阻力問題進行了大量的研究，但關于壁面對顆粒團簇受力影響的研究較少。由于LBM 能實現對流固邊界的精細刻畫及對單個顆粒受力的計算，因此，本文基于LBM 研究顆粒團簇在窄縫壁面影響下的升力、阻力變化規律，為支撐劑輸送研究提供理論依據。

1 LBM 數值模型

1.1 LBM 模型

格子玻爾茲曼模型包含三個要素[15]：流體粒子的離散速度集合、格子演化方程和平衡態分布函數。其中，離散速度集合eα選用三維19 個速度方向的D3Q19 速度模型，如下：

格子演化方程描述了具有離散速度的流體粒子分布函數fi(r,t) 在一個固定格子上不同t時刻的運動過程：

式中：r/m 為格子上節點坐標矢量；ci/(m/s)為速度； Δt/s 為時間步長； Δi為碰撞算子。

碰撞模型采用LADD 等[16]提出的雙松弛模型，對一個立方體的網格來說：

式中：ωi為權重系數，無量綱；cs/(m/s)為理想氣體方程中的聲速；1 是單位矩陣。對D3Q19 模型，ω0=1/3，ω1-6=1/18，ω7-18=1/36，c2s=c2/3。

式中： ρ/(kg/m3)為流體密度；j/(kg/(m2·s))為動量密度；u/(m/s)為速度矢量。

由質量守恒和動量守恒得：

平衡態分布函數常用形式如下[18]：

1.2 邊界條件

固體壁面邊界采用半步長反彈格式[19]，固體邊界格點不在流體網格點上，而是位于流體網格中間，即 (xf+xb)/2處，其表達式為：

LADD[20]對于固液兩相流的移動邊界，提出基于半步長反彈格式加入反映顆粒運動速度的項。假定流體質點位置r恰好位于顆粒表面的外部，并選定離散速度方向cb使r+cbΔt位于顆粒內部格點處，從而對移動邊界進行處理。

流體在流體點和固體點連線的中點發生碰撞，將這些中點依次連接，曲折的邊界線隨著網格的細化逐漸接近顆粒的實際表面邊界線，如圖1 所示。

圖1 不同位置球形顆粒邊界Fig.1 Spherical particle boundary at different positions

1.3 流體-顆粒作用力

采用動量交換法計算顆粒受力[20]，顆粒表面的局部速度可以表示為：

局部速度由顆粒線速度U、角速度 Ω和顆粒質心r確定，其中rb=r+1/2ciΔt是顆粒邊界點的位置。

移動邊界處理導致在流體點和固體點之間產生局部動量交換，但顆粒和流體總的動量保持守恒。邊界上流體對固體顆粒的作用力可以通過動量交換計算出來：

流體作用在顆粒上總的合力與合力矩通過對邊界格點上的f(rb)與rb×f(rb)求和累積得到。

單個球形顆粒的繞流阻力與升力表示為：

式中：CD為曳力系數，無因次；CL為升力系數，無因次； ρ/(kg/m3)為流體的密度；U∞/(m/s)為來流的速度；A/m2為固體在來流方向的投影表面積；D/m 為顆粒直徑；B/m 為壁面縫寬。

1.4 模型求解

采用LADD[20]教授開發的懸浮顆粒流開源軟件Susp3D 進行問題求解。LBM 模擬需要將物理單位轉換到格子單位[21-22]進行求解計算，且顆粒間相對位置固定以保證顆粒團簇形狀參數不變，使用周期邊界條件以有效地模擬顆粒團簇在無限靜止流體中平行于壁面的等速運動。其中：Re=0.1、1 時，計算域為100×27×100；Re=10、100 時，計算域為100×62×100。模擬循環1500 個，共計計算150 000個格子時間步，從而獲取顆粒曳力、升力。

計算采用Intel E5-2670 處理器，工作主頻2.3 GHZ，24 線程并行。若在120 h 的計算時間后，顆粒受力未達到穩定狀態，則通過調整顆粒速度、流體粘度、顆粒粒徑的參數重新計算至收斂。單顆粒壁面影響下(顆粒中心到壁面距離與顆粒粒徑比值H/d=1)曳力收斂曲線如圖2 所示。

圖2 曳力收斂曲線圖Fig.2 Drag convergence curve

2 顆粒團簇物理模型

通過高清相機觀測支撐劑運移過程中顆粒微觀結構[23]如圖3 所示，存在多個顆粒緊密連接聚集成團[24-25]的現象，形成不同構型的支撐劑顆粒團簇。

圖3 顆粒團簇微觀結構Fig.3 Microstructure of particle clusters

根據顆粒團簇運移迎流面的形狀選取三種具有代表性的構型：星型、三棱柱型、長方體型，如圖4 所示。顆粒團簇在無限靜止流體中沿x軸正方向，平行于壁面做等速運動。

圖4 顆粒團簇構型Fig.4 Configuration of particle clusters

3 模擬結果及分析

3.1 模型驗證

不同雷諾數Re下星型顆粒團簇在無壁面流場中運移的曳力系數Cd與BEETSTRA等[5]的模擬結果、TRAN-CONG 等[26]的實驗結果以及CLIFT 等[27]和STOKES[28]的曳力系數公式相吻合如圖5 所示，從而驗證了模型的準確性。對比發現，大雷諾數時模擬結果與文獻存在一定的偏差，這是由于后顆粒受到沿來流方向前顆粒尾流的影響[29]以及顆粒間距的變化會干擾顆粒受力[30]。

圖5 無壁面時團簇曳力系數對比Fig.5 Comparison of drag coefficients of particle cluster without wall

模擬與ZENG 等[31]相同的工況下(顆粒中心到壁面距離與顆粒粒徑比值H/d=1)，有壁面影響時，單顆粒曳力系數Cd與Re關系如圖6 所示，通過與ZENG 等[31]的數值模擬結果對比可以發現，最大誤差小于1.8%，從而驗證了在壁面影響下該模型的有效性。

圖6 壁面影響下單顆粒曳力系數對比Fig.6 Comparison of drag coefficient of single particle with wall

3.2 壁面影響下顆粒團簇曳力與升力模擬

壁面會影響顆粒團簇受力，從而對顆粒沉降和運移產生干擾作用，在狹縫中這種現象更為明顯。為研究此問題，在垂直于y方向增加一組平行壁面，如圖4 所示，并保持縫間距為5 倍顆粒粒徑。模擬壁面影響下，不同Re(顆粒團簇位于壁面中心)、不同顆粒團簇與壁面相對距離H/d時顆粒團簇的運移(H/d是最外側顆粒中心到壁面距離與顆粒粒徑比值，星型及三棱柱型取0.7、0.9、1.1、1.3、1.4，長方體型取0.7、0.9、1.1、1.5、1.9)，從而獲取顆粒的曳力系數、升力系數。

3.2.1 星型團簇

星型團簇曳力系數CD與雷諾數Re的關系如圖7 所示。與單顆粒曳力系數相同，CD隨著Re增大，CD逐漸減小。

圖7 不同雷諾數時星型團簇曳力系數Fig.7 Drag coefficient of star clusters at different Re

不同雷諾數Re時星型顆粒團簇各顆粒曳力系數與團簇平均曳力系數的比值CD/CD,s，如圖8 所示。中心顆粒和前顆粒處于流場中心位置，與單顆粒曳力相似，其曳力主要受來流Re影響，CD隨Re增大而減小，CD/CD,s隨之減小。雷諾數較小時，后顆粒隨Re的增大受壁面邊界影響劇烈，CD/CD,s隨之增大；雷諾數增大至100 時，壁面邊界效應對后顆粒影響減弱，曳力主要受流場Re影響，CD/CD,s隨Re增大而減小。上下左右(近壁面)四顆粒曳力受壁面影響，CD/CD,s隨Re增大而增大。就各顆粒間的曳力大小而言，壁面影響下近壁面顆粒曳力占團簇總曳力比重較大，中心顆粒占比小，說明壁面不僅會加劇對外側顆粒曳力受力的影響，且外側顆粒對內部顆粒的受力有屏障作用。星型團簇各顆粒升力系數CL與雷諾數Re的關系，如圖9 所示。隨著Re的增大，CL總體呈下降趨勢；壁面影響下，邊界效應對上下左右(近壁面)顆粒影響作用明顯，CL數量級較大；由于垂直于z軸方向未設置壁面，上下顆粒無壁面(z)方向CL總小于壁面(y)方向CL；其它顆粒位于流場中心，受到的升力小及對應數量級較小，CL基本保持不變。

圖8 不同雷諾數時星型團簇各顆粒曳力系數與團簇平均曳力系數比值Fig.8 Ratio of each particle drag coefficient of star cluster to average drag coefficient of cluster at different Re

圖9 不同雷諾數時星型團簇各顆粒升力系數Fig.9 Lift coefficient of each particle of star cluster at different Re

星型顆粒團簇不同H/d時曳力系數CD，如圖10所示。隨著H/d的增大，壁面影響逐漸減小，顆粒團簇CD逐漸減小；在H/d=0.7～H/d=0.9 階段，由于壁面邊界層效應的影響較為明顯，顆粒團簇在遠離邊界層的過程中CD下降較快；而在H/d=0.9～H/d=1.4 階段，顆粒團簇從近壁面流場遠離至中心流場處，在這一區域距離邊界層較遠，受壁面影響較小，CD下降程度較小。不同H/d時星型團簇各顆粒曳力系數與團簇平均曳力系數的比值CD/CD,s，如圖11 所示。上下及右顆粒隨H/d增大逐漸靠近壁面，邊界效應增強，CD/CD,s隨之增大；左(近壁面)顆粒隨H/d增大逐漸遠離壁面，邊界效應CD/CD,s隨之減小；H/d=1.4時，左右兩顆粒處于流場對稱位置，CD/CD,s相等；中心顆粒及前后顆粒，曳力主要受來流的流場影響，CD/CD,s隨H/d增大保持平穩且占比較小，說明外側顆粒對內部顆粒受力有屏蔽作用。

圖10 不同H/d 時星型團簇曳力系數Fig.10 Drag coefficient of star clusters at different H/d

圖11 不同H/d 時星型團簇各顆粒曳力系數與團簇平均曳力系數比值Fig.11 Ratio of each particle drag coefficient of star cluster to average drag coefficient of cluster at different H/d

不同H/d時星型顆粒團簇各顆粒壁面方向升力系數CL，如圖12 所示。團簇內除前后顆粒外所有顆粒在周圍顆粒間屏蔽作用下，CL都隨著H/d的變大而小幅度減小；前后側顆粒由于沒有其它顆粒的屏蔽作用，受壁面邊界效應影響明顯，CL大幅度減小。

圖12 不同H/d 時星型團簇各顆粒壁面方向升力系數Fig.12 Lift coefficient of each particle wall direction of star cluster at different H/d

3.2.2 三棱柱型

三棱柱型顆粒團簇曳力系數CD與雷諾數Re關系，如圖13 所示，壁面影響下CD隨Re增大而減小。

圖13 不同雷諾數時三棱柱團簇曳力系數Fig.13 Drag coefficient of triangular prism cluster under different Re

不同雷諾數Re時三棱柱顆粒團簇各顆粒曳力系數與團簇平均曳力系數比值CD/CD,s，如圖14 所示。中心及前顆粒處于流場中心位置，壁面影響較弱，曳力主要與Re有關，CD/CD,s隨Re增大而減小；后側近壁面顆粒在邊界效應影響下，CD/CD,s隨Re增大而增大；后排中心顆粒，在小雷諾數時CD受壁面影響先增大，Re較大時CD主要受流場Re影響驟減，CD/CD,s隨之先增大再減小；就各顆粒曳力大小而言，由于壁面影響下外側顆粒對內部顆粒受力有屏障作用，后側近壁面顆粒曳力比重較其他顆粒大。

圖14 不同雷諾數時的三棱柱團簇各顆粒曳力系數與團簇平均曳力系數比值Fig.14 Ratio of each particle drag coefficient of triangular prism cluster to average drag coefficient of cluster at different Re

三棱柱型顆粒團簇升力系數CL與雷諾數Re的關系，如圖15 所示，前顆粒及后排中間顆粒位于流場中心，受到升力及相應CL量數級較小。其它顆粒位于團簇外側，壁面(y)方向及無壁面(z)方向升力系數CL均隨Re的增大而減小；且壁面方向CL總小于無壁面方向CL。

圖15 不同雷諾數時三棱柱團簇升力系數Fig.15 Lift coefficient of triangular prism cluster under different Re

三棱柱型顆粒團簇不同H/d下曳力系數CD與雷諾數Re的關系，如圖16 所示。壁面影響下三棱柱型顆粒團簇CD隨著H/d的增大而減小；在H/d=0.7～H/d=0.9 階段，顆粒團簇逐漸遠離壁面，由于邊界層效應的影響較為明顯，CD下降較快；在H/d=0.9～H/d=1.4 階段，顆粒團簇遠離壁面至中心流場處，距邊界層較遠，受壁面影響較小，CD下降程度較小。

圖16 不同H/d 時三棱柱團簇曳力系數Fig.16 Drag coefficient of triangular prism cluster at different H/d

三棱柱型顆粒團簇不同H/d時，各顆粒曳力系數與團簇平均曳力系數比值CD/CD,s，如圖17 所示。隨H/d的增大，后排近壁面顆粒逐漸遠離壁面，邊界效應減弱，CD/CD,s顯著減小；后排遠壁面顆粒逐漸靠近壁面，邊界效應加強，CD/CD,s顯著增大；H/d=1.4 時，顆粒團簇位于流場中間位置，后排近壁面與遠離壁面顆粒處于流場對稱位置，CD/CD,s相同；在外側顆粒的屏蔽作用下，中間及后排中間顆粒CD/CD,s保持平穩；壁面邊界效應會加劇顆粒受力，使后排近壁面顆粒曳力占團簇總曳力比重較大。三棱柱型顆粒團簇不同H/d下各顆粒壁面方向升力系數CL結果，如圖18 所示。分析發現，不同H/d時，近壁面顆粒隨H/d的增大逐漸遠離邊界層壁面方向CL減小，遠壁面顆粒隨著H/d的增大逐漸靠近壁面邊界，壁面方向CL增大；壁面邊界效應加劇顆粒受力，使近壁面顆粒CL總大于遠壁面顆粒CL。

圖17 不同H/d 時三棱柱團簇各顆粒曳力系數與團簇平均曳力系數比值Fig.17 Ratio of each particle drag coefficient of triangular prism cluster to average drag coefficient of cluster at different H/d

圖18 不同H/d 時三棱柱團簇各顆粒壁面方向升力系數Fig.18 Lift coefficient of each particle wall direction of triangular prism cluster at different H/d

3.2.3 長方體型

長方體型顆粒團簇曳力系數CD與雷諾數Re的關系，如圖19 所示，與星型及三棱柱型顆粒團簇類似，壁面影響下長方體型顆粒團簇CD隨Re的增大。

圖19 不同雷諾數時長方體型團簇曳力系數Fig.19 Drag coefficient of cuboid clusters under different Re

不同雷諾數Re時長方體型顆粒團簇各顆粒曳力系數與團簇平均曳力系數比值CD/CD,s，如圖20所示。可以發現，前排顆粒曳力受來流Re影響，CD/CD,s隨著Re的增大而減小；中間顆粒CD/CD,s基本保持不變；后排顆粒在Re較小時曳力受壁面邊界效應影響，CD/CD,s先增大，Re較大時對曳力的影響增強，CD/CD,s減小。就各顆粒間的曳力大小而言，中間顆粒受力在外側顆粒屏蔽作用下，曳力系數占比較小。

圖20 不同雷諾數時長方體團簇各顆粒曳力系數與團簇平均曳力系數比值Fig.20 Ratio of each particle drag coefficient of cuboid cluster to average drag coefficient of cluster at different Re

不同雷諾數Re時長方體型顆粒團簇各顆粒升力系數CL，如圖21 所示。通過分析發現，壁面影響下各顆粒壁面(y)方向和無壁面(z)方向CL均隨著Re的增大而減小；且壁面影響下，顆粒的無壁面方向CL總體上小于壁面方向的結果。

圖21 不同雷諾數時長方體型團簇各顆粒升力系數Fig.21 Lift coefficient of each particle of cuboid cluster at different Re

長方體型顆粒團簇不同H/d下曳力系數CD的結果，如圖22 所示。隨著H/d的增大，壁面影響減小，顆粒團簇CD逐漸減小；顆粒團簇在遠離邊界層的H/d=0.7～H/d=0.9 階段，由于壁面邊界層效應的影響較為明顯，CD下降較快；在H/d=0.9～H/d=1.9階段，顆粒團簇從靠近壁面流場處遠離至中心流場，距離壁面邊界層較遠，受壁面影響較小，CD下降程度較小。

圖22 不同H/d 時長方體團簇曳力系數Fig.22 Drag coefficient of cuboid clusters at different H/d

不同H/d時長方體型顆粒團簇各顆粒曳力系數與團簇平均曳力系數比值CD/CD,s，如圖23 所示。隨著H/d的增大，近壁面前中后三個顆粒逐漸遠離壁面，壁面邊界效應減弱，CD/CD,s逐漸減小；遠離壁面前中后三顆粒逐漸靠近壁面，壁面邊界效益增強，CD/CD,s隨之增大；就各顆粒曳力大小而言，邊界效應會加劇近壁面顆粒受力，使遠壁面顆粒CD/CD,s較小；內側顆粒受力在團簇外側顆粒的屏蔽作用下，中間顆粒曳力系數占比較小。

圖23 不同H/d 時長方體型團簇各顆粒曳力系數與團簇平均曳力系數比值Fig.23 Ratio of each particle drag coefficient of cuboid cluster to average drag coefficient of cluster at different H/d

長方體型顆粒團簇不同H/d下各顆粒升力系數CL結果，如圖24 所示。團簇內各顆粒壁面方向CL均隨H/d的變大而減小；近壁面顆粒的CL總大與遠壁面顆粒；中間兩顆粒壁面方向升力在外側顆粒屏蔽作用下，隨著H/d的變大基本保持不變。

圖24 不同H/d 時長方體型團簇各顆粒壁面方向升力系數Fig.24 Lift coefficient of each particle wall direction of cuboid cluster at different H/d

4 結論

本文構建了星型、三棱柱型、長方體型顆粒團簇模型，并基于LBM 研究了窄縫中顆粒團簇在不同雷諾數Re、不同顆粒團簇與壁面相對距離H/d下的升阻力系數變化規律，得到以下結論：

(1)星型、三棱柱型、長方體型顆粒團簇的曳力系數CD隨著Re、H/d的增大而減小；且顆粒團簇在H/d=0.7～0.9 階段，距邊界層較近，受壁面影響明顯，CD下降較快；在顆粒團簇繼續逐漸遠離至中心流場階段，距邊界層較遠，受壁面影響程度小，CD下降較慢。

(2)星型、三棱柱型、長方體型顆粒團簇，顆粒曳力系數占總團簇平均曳力系數比值CD/CD,s，隨Re增大而減小及隨壁面邊界效應減弱而增大的程度不同，前排顆粒曳力系數比值減小、中間顆粒先增大后減小、近壁面顆粒逐漸增大。隨著H/d的增大，近壁面顆粒遠離壁面使邊界效應減弱，比值隨之減小；遠壁面顆粒反之；中間顆粒位于流場中心，比值保持穩定。

(3)星型、三棱柱型、長方體型顆粒團簇，各顆粒升力系數CL隨著Re、H/d的增大總體呈減小趨勢；且壁面影響下隨Re的增大，顆粒壁面方向CL均小于無壁面方向的結果。

(4)星型、三棱柱型、長方體型顆粒團簇，近壁面顆粒曳力系數與總團簇平均曳力系數比值CD/CD,s及升力系數CL，總大于遠壁面顆粒的結果；且外側顆粒受力總大于中間顆粒。壁面影響下，邊界效應會加劇近壁面顆粒受力，且顆粒團簇外側顆粒對內部顆粒受力有屏蔽作用。