基于XDLVO理論的黏性泥沙絮凝模擬格子玻耳茲曼模型

喬光全，張慶河，張金鳳，程洪劍，盧 昭

(1. 天津大學(xué)水利工程仿真與安全國家重點實驗室，天津 300072；2. 中交水運規(guī)劃設(shè)計院有限公司，北京 100007)

黏性泥沙在我國海岸河口地區(qū)分布廣泛，其絮凝沉降過程與港口和航道淤積、海岸河口環(huán)境變化等密切相關(guān)，對黏性泥沙絮凝沉降特性的研究具有十分重要的意義．海岸河口地區(qū)影響?zhàn)ば阅嗌承跄两档囊蛩丶扔心嗌愁w粒礦物質(zhì)組成、泥沙濃度等黏性泥沙自身特性，又有海水鹽度和河口淡水屬性等水體環(huán)境特性．因此，描述黏性泥沙絮凝沉降時，需要盡可能包括上述各種因素的影響．

近年來，一些學(xué)者開始利用數(shù)值模擬方法研究黏性泥沙絮凝沉降，如洪國軍等[1]用受限的絮團聚集模型(diffusion limited cluster aggregation，DLCA)研究了黏性細(xì)顆粒泥沙的絮凝沉降過程，得到了分形維數(shù)和顆粒濃度的關(guān)系，討論了絮凝-沉降過程的各種特性．王龍等[2]將 XDLVO 理論(extended Derjaguin-Landau-Verwey-Overbeek theory)和無參數(shù)可解性的加速斯托克斯動力學(xué)(accelerated Stokesian dynamics)方法結(jié)合，對黏性泥沙絮凝沉降進(jìn)行了數(shù)值研究，獲得了鹽度、濃度和水合作用參數(shù)對泥沙沉降絮凝的影響．上述數(shù)值模型雖然能夠獲得有關(guān)黏性泥沙絮凝變化的規(guī)律，但不能完整描述顆粒絮凝的水動力過程，因而難以全面考慮各種因素，如水流紊動等對絮凝的影響．

格子玻耳茲曼(lattice Boltzmann，LB)方法是一種近年來得到快速發(fā)展的流動模擬方法，具有編程簡單、易于并行求解顆粒與流體全尺度相互作用問題等優(yōu)點．在合理考慮顆粒之間相互作用力和定義顆粒絮凝狀態(tài)的情況下，LB方法可以從微觀水動力角度完整描述絮凝過程，為深入研究黏性泥沙絮凝規(guī)律提供了工具．張金鳳等[3-4]采用LB方法，考慮顆粒之間的范德華吸引力，建立了黏性泥沙不等速沉降絮凝的數(shù)值模型．盧昭[5]和程洪劍[6]在該模型中引入了DLVO(Derjaguin-Landau-Verwey-Overbeek)理論，考慮了顆粒間范德華力和雙電層作用力，模擬了鹽度和濃度對絮凝的影響．上述工作逐步合理地考慮了顆粒間相互作用力，但由于模型中沒有反映泥沙礦物種類的參數(shù)，因此尚未能對礦物種類等因素對黏性泥沙絮凝的影響進(jìn)行深入分析．

為了在 LB模型中更全面地反映各種因素對黏性泥沙絮凝沉降的影響，筆者將顆粒之間相互作用的XDLVO理論引入 LB模型，并利用模型分別對伊利土、蒙脫土和高嶺土3種黏土礦物泥沙顆粒間的絮凝進(jìn)行模擬，以分析礦物質(zhì)成分對絮凝的影響．

1 黏性泥沙絮凝的LB模型

1.1 LB模型簡介

LB方法采用微觀流體的粒子速度分布函數(shù)描述流體運動，fi的演化方程為

式中：fi為在離散時間t內(nèi)、格子r處速度為 ci的假想流體粒子數(shù)；i為速度方向；Δi為線性碰撞算子，表示為

LB模型可模擬顆粒和流體的相互作用．Ladd[7]提出了處理顆粒運動的LB算法：流體格點位于顆粒的外部，每個格點的速度為bc，bt+Δrc位于顆粒的內(nèi)部，格子的中點表示固體顆粒邊界，格點的連線穿過流-固邊界表面，流體粒子沿著格子連線運動并與懸浮固體顆粒發(fā)生碰撞和相互作用，在流體和固體間的局部能量交換過程中保持能量守恒．流體粒子每個時間步離開原來格點沿規(guī)定方向運動到鄰近格點，當(dāng)流體粒子運動到固定壁面時，利用全反射條件可滿足壁面無滑移．因此，流體-固體顆粒邊界的離散模型為

顆粒的平動速度和轉(zhuǎn)動角速度可根據(jù)平動加速度α和轉(zhuǎn)動角加速度 β表達(dá)式逐步積分獲得，其表達(dá)式分別為

式中：m和 J分別為泥沙顆?；蛐鯃F的質(zhì)量和動量矩；F和M分別為泥沙顆?；蛐鯃F受到的總外力和總力矩．對于黏性細(xì)顆粒泥沙，顆粒受到的總力可表示為

1.2 顆粒間相互作用力

20世紀(jì)80年代以來，Somasundaran[9]和van Oss等[10]在研究黏土礦和磁鐵礦的凝聚行為時發(fā)現(xiàn)，極性溶液中距離很近的 2個顆粒表面由于氫離子的黏結(jié)作用會形成路易斯酸堿能(acid-base energy)，且在短距離情況下，這種 A-B相互作用能較靜電排斥能和范德華吸引能大 1～2個量級．將這種作用加入到DLVO理論中，形成了 XDLVO理論．胡岳華等[11]也認(rèn)識到這種力的存在，給出了顆粒間相互作用能的計算公式，并通過物理實驗驗證了XDLVO理論的適用性．Wu等[12]用XDLVO理論和DLVO理論分別對各種礦物顆粒懸浮液的穩(wěn)定性進(jìn)行了研究，發(fā)現(xiàn)前者比后者得到的結(jié)果更為精確．

由于 XDLVO理論考慮了顆粒之間的范德華吸引力、靜電排斥力和A-B作用力，更為合理地描述了細(xì)顆粒及膠體在電解質(zhì)溶液中的凝聚行為，因此這里采用 XDLVO理論描述黏性泥沙顆粒間的相互作用．XDLVO理論將顆粒之間的總勢能表示為

1.2.2 雙電層作用力

式中：e為元電荷，e＝1.6×10-19,C；NA為阿伏伽德羅常數(shù)，NA＝6.022×1023/mol；k為 Boltzmann 常數(shù)，k＝1.38×10-23,J/K；T為絕對溫度，K；cm為陽離子濃度，mol/L；zm表示陽離子化合價，無量綱．

雙電層作用力為雙電層作用能對距離的負(fù)導(dǎo)數(shù)，即

1.2.3 A-B作用力

A-B作用能是XDLVO理論區(qū)別于DLVO理論的重要作用能，其表達(dá)式[14]為

式中：λAB為衰減長度，與顆粒本身的特性有關(guān)，具有

1.3 LB模型中XDLVO勢能合理性分析

A-B力中較為敏感的參數(shù)為衰減長度ABλ[3]，為了驗證模型的合理性，這里模擬不同衰減長度條件下，2個不等直徑的球形泥沙顆粒在靜水中的沉降過程．計算范圍取為 0.1,mm×0.2,mm×0.1,mm，格子數(shù)為 50×100×50，2個泥沙顆粒直徑分別為 10,μm和 5,μm．在計算中保持顆粒不絮凝，并使之盡量靠近．圖 1顯示了衰減長度取不同值時顆粒間各種勢能，圖 1中的曲線表示公式計算值，離散點表示 LB模型模擬結(jié)果．從圖中可以看出，模擬結(jié)果和公式計算值完全吻合．隨著ABλ增大，A-B勢能的作用范圍與數(shù)值也增大，XDLVO總勢能的勢阱深度減小甚至消失，勢壘高度增加，同時總勢能的勢阱位置和勢壘位置都略有變大．該結(jié)果與文獻(xiàn)[3]和文獻(xiàn)[11]的模擬分析結(jié)果一致，說明LB模型中引進(jìn)的XDLVO勢是合理的．

圖1 顆粒間勢能和衰減長度的關(guān)系Fig.1 Relationship between the potential and decay length

2 黏土礦物對黏性泥沙絮凝影響的模擬分析

2.1 模擬參數(shù)設(shè)置

選取伊利土、蒙脫土和高嶺土3種常見黏土礦物進(jìn)行模擬分析，以檢驗基于 XDLVO理論所建立的LB模型是否反映了不同黏土礦物的絮凝特性．3種物質(zhì)和水的相關(guān)參數(shù)見表 1[13-15]，計算中，水和顆粒密度分別為 1,000,kg/m3和 2,650,kg/m3，水的黏滯系數(shù)為 1.0×10-6,m2/s，重力加速度為 9.8,m/s2．模擬算例考慮重力作,用下做靜水沉降的 2個泥沙顆粒，直徑分別為10,μm和5,μm，初始時刻大顆粒位于小顆粒上方，兩者垂向間距為2,μm，水平方向中心位置偏心距在0～4,μm間變動．

表1 計算參數(shù)Tab.1 Parameters in simulation

2.2 模擬結(jié)果分析

當(dāng)兩泥沙顆粒表面帶有的束縛水層(滑動層)相接觸，形成公共滑動層時，可認(rèn)為顆粒結(jié)合而絮凝[16]，因此計算時以 2倍滑動層厚度作為泥沙絮凝的判斷標(biāo)準(zhǔn)．圖 2顯示了不同算例中絮團的形成時間，圖 3顯示了沉降過程中顆粒的運動軌跡．由圖2和圖3可知：初始位置相同時，不同泥沙礦物發(fā)生絮凝的時間不同，伊利土最先發(fā)生絮凝，其次是高嶺土和蒙脫土；泥沙礦物組成相同時，顆粒的初始位置影響其絮凝過程，兩泥沙顆粒中心無偏離時，絮凝發(fā)生最早，偏心距離增大，絮凝出現(xiàn)時間延后，高嶺土和蒙脫土在泥沙顆粒偏心距離為 3,μm 時不再發(fā)生絮凝，大顆粒繞過小顆粒后分別獨自沉降，伊利土達(dá)到不絮凝的顆粒初始偏心距離為 4,μm．由此可知，就本文所取的3種黏土礦物而言，絮凝程度從易到難的順序依次是伊利土、高嶺土和蒙脫土．

圖2 不同礦物成分的泥沙顆粒絮凝時間Fig.2 Flocculation time of sediment particles for different minerals

圖3 顆粒運動路徑Fig.3 Trajectories of particle settling

結(jié)合XDLVO理論，對沉降過程中泥沙顆粒的運動軌跡和受力進(jìn)行分析，可以對上述過程做出合理的解釋．以泥沙顆粒水平偏心距為 2,μm的情況為例進(jìn)行分析，圖4顯示了3組黏土礦物泥沙顆粒間距的歷時曲線，圖5為圖4的局部放大．圖6顯示了2個泥沙顆粒在沉降過程中的受力歷時曲線，圖中從上至下分別為水平和垂直方向作用力以及顆粒之間的XDLVO力．

兩顆粒在距離較遠(yuǎn)(0.4,s之前)時，主要受到重力和周圍流體運動引起的外力，XDLVO力的作用不明顯，不同黏土顆粒其運動軌跡和受力狀態(tài)基本相同．x方向上，大小顆粒之間表現(xiàn)為相互排斥，且排斥力隨沉降時間增加而不斷增大；y方向上，大顆粒受下方小顆粒尾流的影響，受力更為復(fù)雜，但大小顆粒之間仍表現(xiàn)為相互排斥．

圖4 泥沙顆粒間距離歷時曲線Fig.4 Time series of distance between particles

圖5 泥沙顆粒間距離歷時曲線(局部放大)Fig.5 Time series of distance between particles with local zooming

顆粒沉降 0.4,s后，隨著顆粒距離減小，XDLVO力的作用開始明顯，兩顆粒之間的排斥力迅速減小直至變?yōu)橹饾u增大的吸引力，顆粒間距離迅速變小，3種黏土顆粒先后發(fā)生絮凝．由表 1可以得知，伊利土的A-B自由能為負(fù)值，A-B力表現(xiàn)為吸引力，且其衰減長度大于其他2種礦物質(zhì)，所以XDLVO力作用范圍大，兩個顆粒距離較遠(yuǎn)時顆粒間總作用力就表現(xiàn)為吸引力(見圖5和圖6(a))，因此容易發(fā)生絮凝．對于高嶺土和蒙脫土，A-B力表現(xiàn)為斥力，且作用范圍小，使得顆粒間總的作用力(吸引力)小于伊利土顆粒間作用力，且和伊利土的該力相比，作用距離小，作用時間晚(見圖5、圖6(b)和圖6(c))，從而導(dǎo)致絮凝能力弱于伊利土．高嶺土的 Hamaker常數(shù)大于蒙脫土，其顆粒間范德華吸引力大于蒙脫土；高嶺土的AB自由能參數(shù)小于蒙脫土，其顆粒間 A-B斥力小于蒙脫土．因此，高嶺土顆粒間 XDLVO力影響范圍大于蒙脫土，絮凝時間較早，絮凝能力稍強．

對于同一種黏土礦物，泥沙顆粒初始位置無水平偏移時，兩顆粒對心正碰，容易發(fā)生絮凝．若初始位置水平偏心不大，運動過程中顆粒間最小距離在XDLVO力的作用范圍內(nèi)，且吸引力大于排斥力，合力的作用仍會使泥沙顆粒形成公共滑動層，從而發(fā)生絮凝現(xiàn)象．當(dāng)偏心距離大到兩顆粒運動過程中其最小距離大于XDLVO力的作用范圍，或者吸引力小于排斥力時絮凝現(xiàn)象就不發(fā)生．由于伊利土的 XDLVO力的數(shù)值和作用范圍都較高嶺土和蒙脫土大，因此其不發(fā)生絮凝的偏心距離大于后兩者．

總的來看，圖3～圖6所輸出的模擬結(jié)果充分反映了黏土顆粒在水體中沉降發(fā)生絮凝的變化過程．文獻(xiàn)[17-18]分別通過物理實驗得出伊利土絮凝能力高于高嶺土和高嶺土絮凝能力高于蒙脫土的結(jié)論，并解釋了其絮凝機制．本文的模擬結(jié)果和上述實驗結(jié)果一致，可以認(rèn)為本文所建立的模型能夠從顆粒之間相互作用的微觀動力角度描述顆粒之間的絮凝過程，較好地反映了不同黏土礦物的絮凝特性．

圖6 不同礦物組分的泥沙顆粒受力歷時曲線Fig.6 Time series of forces between particles for different minerals

3 結(jié) 論

(1)考慮XDLVO力的格子Boltzmann模型能夠較好地模擬黏性泥沙的絮凝沉降以及該過程中顆粒之間的受力情況，從微觀水動力學(xué)角度反映了絮凝機理和規(guī)律，可以用于進(jìn)一步研究各種因素對黏性泥沙絮凝的影響規(guī)律．

(2)不同的黏土礦物顆粒絮凝的難易程度不同，就本文模擬的3種黏土礦物而言，絮凝程度從易到難的順序依次是伊利土、高嶺土和蒙脫土，和以往的認(rèn)識一致．

［1］ 洪國軍，楊鐵笙. 黏性細(xì)顆粒泥沙絮凝及沉降的三維模擬[J]. 水利學(xué)報，2006，37(2)：172-177.

Hong Guojun，Yang Tiesheng. 3-D simulation of flocculation-settling of cohesive fine sediment [J].Journal of Hydraulic Engineering，2006，37(2)：172-177(in Chinese).

［2］ 王 龍，李家春，周濟福. 黏性泥沙絮凝沉降的數(shù)值研究[J]. 物理學(xué)報，2010，59(5)：3315-3323.

Wang Long，Li Jiachun，Zhou Jifu. Numerical study of flocculation settling of cohesive sediment [J].Acta Physica Sinica，2010，59(5)：3315-3323(in Chinese).

［3］ 張金鳳，張慶河. 黏性泥沙不等速沉降絮凝的格子Boltzmann模擬[J]. 水利學(xué)報，2009，40(4)：385-390.Zhang Jinfeng，Zhang Qinghe. Lattice Boltzmann simu-lation for flocculation of cohesive sediment due to differential settling [J]. Journal of Hydraulic Engineering，2009，40(4)：385-390(in Chinese).

［4］ Zhang Jinfeng，Zhang Qinghe. Lattice Boltzmann simulation of the flocculation process of cohesive sediment due to differential settling [J].Continental Shelf Research，2011，31(10)：94-105.

［5］ 盧 昭. 鹽度對黏性泥沙絮凝影響的格子波耳茲曼模擬研究[D]. 天津：天津大學(xué)建筑工程學(xué)院，2007.

Lu Zhao. Lattice Boltzmann Simulation for Effects of Salinity on Flocculation of Cohesive Sediment [D].Tianjin：School of Civil Engineering，Tianjin University，2007(in Chinese).

［6］ 程洪劍. 黏性泥沙不等速沉降絮凝的格子波耳茲曼模擬研究[D]. 天津：天津大學(xué)建筑工程學(xué)院，2008.

Cheng Hongjian. Lattice Boltzmann Simulation for Differential Settling Flocculation of Cohesive Sediment[D]. Tianjin：School of Civil Engineering，Tianjin University，2008(in Chinese).

［7］ Ladd A J C. Numerical simulations of particulate suspensions via a discretized Boltzmann equation(Part 1)：Theoretical foundation [J].Journal of Fluid Mechanics，1994，271：285-309.

［8］ Nguyen N Q，Ladd A J C. Lubrication corrections for lattice-Boltzmann simulations of particle suspensions[J].Physical Review E-Statistical， Nonlinear，and Soft Matter Physics，2002，66(4)：1-12.

［9］ Somasundaran P. Role of surface phenomena in the beneficiation of fine particles[J]. Mining Engineering，1984，36(8)：1177-1186.

［10］ van Oss C J，Chaudhury M K，Good R J. The mechanism of phase separation of polymers in organic mediaapolar and polar systems[J].Separation Science and Technology，1989，24(1/2)：15-30.

［11］ 胡岳華，邱冠周，王淀佐. 細(xì)粒浮選體中擴展的DLVO理論及應(yīng)用[J]. 中南礦冶學(xué)院學(xué)報，1994，25(3)：310-314.

Hu Yuehua，Qiu Guanzhou，Wang Dianzuo． Extended DLVO theory and its applications in flotation of fine particles[J].Journal of Central South Institute of Mining and Metallurgy，1994，25(3)：310-314(in Chine-se).

［12］ Wu W，Giese R F，van Oss C J. Stability versus flocculation of particle suspensions in water-correlation with the extended DLVO approach for aqueous systems，compared with classical DLVO theory[J].Colloids and Surfaces B：Biointerfaces，1999，14(1/2/3/4)：47-55.

［13］ Hoek E M V，Agrwal G K. Extended DLVO interactions between spherical particles and rough surfaces[J].Journal of Colloid and Interface Science，2006，298(1)：50-58.

［14］ Li Z，Giese R，van Oss C J. Surface thermodynamic properties of some illites[C]//Program and Abstracts for Clay Minerals Society，28th Annual Meeting. Houston，TX：the Lunar and Planetary Institute，1991：100.

［15］ Duran J D G，Ramos-Tejada M M，Arroyo F J，et al.Rheological and electrokinetic properties of sodium montmorillonite suspensions(I)：Rheological properties and interparticle energy of interaction[J]. Journal of Colloid and Interface Science，2000，229(1)：107-117.

［16］ 楊鐵笙，熊祥忠，詹秀玲，等. 黏性泥沙懸浮液中顆粒表面滑動層厚度的計算[J]. 水利學(xué)報，2002(5)：20-25.

Yang Tiesheng，Xiong Xiangzhong，Zhan Xiuling，et al. The study on slipping water layers of cohesive sediment particles [J].Journal of Hydraulic Engineering，2002(5)：20-25(in Chinese).

［17］ 王毓華，陳興華，胡業(yè)民，等. 磷酸鹽對細(xì)粒鋁硅酸鹽礦物分散行為的影響[J]. 中南大學(xué)學(xué)報：自然科學(xué)版，2007，38(2)：238-244.

Wang Yuhua，Chen Xinghua，Hu Yemin，et al. Influences of phosphates on dispersion of fine aluminumsilicate minerals [J].Journal of Central South University：Science and Technology，2007，38(2)：238-244(in Chinese).

［18］ Tombacez E，Szekeres M. Surface charge heterogeneity of kaolinite in aqueous suspension in comparison with montmorillonite [J].Applied Clay Science，2006，34(1/2/3/4)：105-124.