分形結構團聚體在流場中分散混合的數值模擬*

麻向軍 邱忠財 陶成文

(華南理工大學聚合物新型成型裝備國家工程研究中心，廣東廣州510640)

分散混合在材料加工中起著重要作用，尤其是顆粒團聚體在聚合物熔體中的混合與復合過程．分散混合的目的是在流體拖曳力作用下，打破分散相顆粒團聚體內的物理鍵合，使其發生剝蝕或破碎，從而獲得更小的碎片或者顆粒［1］．預測和控制分散混合過程對設計和操作混合設備進而提高最終制品性能有重要意義，許多研究者從理論和實驗兩方面對分散混合現象進行了研究．Feke等［2］提出了一種用于確定團聚體在剪切流場中破碎的概率模型;Horwatt等［3］將處于剪切流場的非均質團聚體視為一個內聚力作用下的連續體，確定了團聚體破碎時的臨界流體拖曳力;Scurati等［4］考慮了團聚體尺寸、流體黏度和剪切速率對分散過程的影響，開發了一種用于預測分形結構團聚體在穩態和非穩態剪切流場中分散的剝蝕動力學模型，并與實驗結果進行了對比．

隨著計算機技術的發展，數值模擬廣泛用于研究分散混合現象，常用的模型有連續介質模型和離散介質模型．由于顆粒團聚體是典型的離散體，力學特征“散”而“動”(“散”是指顆粒性質、尺寸和形貌的分散性，“動”是指顆粒運動的瞬時性和非線性)，不滿足連續模型的小變形和位移協調的假設，因此，應采用離散介質模型研究分散混合過程［5］．目前，離散介質模型中應用最廣泛的是Cundall［6］提出的離散單元法(DEM)，其基本原理是將離散介質視為獨立的塊體或顆粒單元的集合，單元之間的相互作用力可根據力和位移的關系確定，單元的運動由牛頓運動定律確定，運動過程中，單元之間可接觸也可分開，允許單元發生很大的平動和轉動，可以求解一些離散介質的力學行為，如堆積、流動和分離［7-9］．近年來，國外一些學者將DEM應用于研究分散混合現象，顯示出這一方法的優越性［10-12］，但是對于分形結構團聚體分散混合的研究較少．文中在經典離散元理論模型中引入范德華力和流體拖曳力，編制了二維離散元程序DEMix2D，并對二維情況下由均質的納米級球形顆粒組成的具有分形結構的團聚體在純剪切、純拉伸和剪切-拉伸混合流場中的分散混合行為進行數值模擬．

1 數值模型

1．1 顆粒運動控制方程

由于分散混合過程非常復雜，必須做一些簡化和假設:①忽略顆粒重力和布朗運動;②顆粒密度和流體密度相近，忽略浮力的作用;③流體處于層流狀態，且只考慮流體對顆粒的作用力，不考慮顆粒對流體的作用力，并且認為顆粒在流體中的任意分布不會影響流場分布;④流場無邊界．滿足以上假設后，團聚體在分散過程中，每個顆粒受到的力包括顆粒間的接觸力、范德華力以及流體拖曳力．顆粒受到這些力作用后的運動可由牛頓第二定律求得:

式中，m、I、vi和ωi分別為顆粒i的質量、轉動慣量、速度和角速度，Fi和Mi分別為顆粒i受到的合力和合力矩．Fi和Mi可表示成:

1．2 接觸模型

顆粒間的接觸力和阻尼力可以根據力與位移的關系(即接觸模型)求得．文中采用Cundall的線性粘彈性模型［10］，即顆粒間可以理解為線性彈簧、線性粘壺和滑片相連接的形式，如圖1所示．

線性粘彈性接觸模型中，接觸力和阻尼力表示為

式中:kn和kt分別為彈簧的法向剛度和切向剛度，kn= RE，kt=kn/10，其中R為顆粒半徑，E為顆粒的彈性模量［10］;dn和 dt分別為法向和切向阻尼系數分別為法向和切向相對位移，vn，ij和 vt，ij分別為法向和切向相對速度．同時，和還須滿足庫倫-摩爾摩擦準則［13］．

1．3 范德華力

根據Hamaker理論，兩個半徑相同的球形顆粒i、j之間的范德華力fijv為

式中，A和h分別表示Hamaker常數和兩顆粒間隔距離，nij為單位向量．

式中，xi和xj為顆粒i和j的位置矢量．

1．4 流體拖曳力

流體拖曳力的計算是模擬顆粒團聚體分散行為的關鍵，目前普遍采用的方法是根據斯托克斯定律求解拖曳力和力矩．顆粒在剪切流場中受到的流體拖曳力如圖2所示．

圖1 顆粒間的接觸模型Fig．1 Contact model between particles

圖2 顆粒在剪切流場中的受力示意圖Fig．2 Schematic diagram of forces on particles in shear flow field

根據斯托克斯定律，單個顆粒處于流場時受到的拖曳力和力矩分別為

式中:η表示流體黏度;e和ω0分別表示流場的速度梯度張量和渦度矢量，與流場類型有關．

對于二維純剪切流場:

式中，γs為剪切速率．

對于單軸拉伸流場:

式中，γe為拉伸速率．

對于剪切-拉伸混合流場，e和ω0可認為是兩種流場的疊加［15］:

兩種流場疊加時形變速率張量第二不變量I2相等，因此，γs和 γe必須滿足 γs=

團聚體中的顆粒與單個顆粒受到的流體拖曳力不同．團聚體中的顆粒由于受到周圍顆粒的影響，計算團聚體中每個顆粒受到的流體拖曳力和力矩時需要引入一個有效因子ζ以反映周圍顆粒的影響．團聚體中顆粒半徑一致時，一個顆粒周圍最多有6個接觸顆粒，此時，該顆粒被接觸顆粒包圍，可認為它沒有受到流體的作用，即ζ=0;當這個顆粒周圍有5個接觸顆粒時，則ζ=1/6;依此類推，當這個顆粒沒有接觸顆粒時，ζ=1．

1．5 運動方程求解

采用基于軟球模型的離散元方法，對運動方程采用歐拉格式求解［16］．任一時刻t，根據團聚體中任一顆粒i在空間的位置及其周圍顆粒的位置，可由接觸模型、范德華力模型和流體拖曳力模型計算出該顆粒受到的接觸力、阻尼力、范德華力和流體拖曳力之和Fi(xi(t))，然后根據運動方程可計算得到t+Δt時刻顆粒i的位置和速度，即

當初始條件給定時，就可以由前一時刻顆粒i的位置和速度確定下一時刻的位置和速度，進而確定團聚體在流場作用下的演變過程．

由于對運動方程采用歐拉格式求解，因此是條件穩定的，時間步長必須小于臨界值才能保證算法穩定．文中采用簡諧振動法確定臨界時間步長，即

式中，C為時間步長安全系數，k為彈簧的剛度．

2 模擬條件及參數

2．1 顆粒團聚體結構生成

在已有的研究中，通常以致密堆積的團聚體為對象，而實際的團聚體結構非常復雜，通常具有分形結構．在二維情況下，以團聚體中心的顆粒為圓心，任意長度r為半徑的圓內含有的顆粒數為N(r)，若N(r)與r滿足以下關系:

則認為團聚體滿足分形結構，其中z為分形維數，它反映了復雜形體占有空間的有效性．根據這個規定，生成團聚體的方法為:首先在流場中心生成一個初始顆粒1，然后生成任意一個與顆粒1相切的顆粒2，再生成第3個顆粒，該顆粒與之前生成的至少1個顆粒相切且不與其他顆粒交叉，同時滿足式(18)，如此重復下去直到生成所需要的顆粒數．圖3為采用上述方法生成的團聚體，顆粒數為130，z為1．8．

2．2 計算參數

模擬過程中，流場強度必須足夠大才能打破團聚體的物理鍵合使其分散．定義流場強度:

圖3 團聚體結構Fig．3 Structure of agglomerates

選取不同流場強度下的剪切、拉伸和剪切-拉伸混合流場進行團聚體分散行為的模擬計算，流場強度分別為460、820和1200Pa．模擬計算中的其他參數如表1所示．

表1 模擬計算中的參數Table 1 Parameters used in simulation

3 模擬結果與分析

3．1 流場形式對分散過程的影響

流場強度為820Pa時，模擬計算得到的團聚體在剪切流場、拉伸流場、剪切-拉伸混合流場中的分散過程如圖4所示．

由圖4可見，3種形式的流場中，團聚體均先發生相應的變形，然后再分散成碎片或者顆粒，團聚體在分散過程中存在剝蝕和破裂兩種分散方式，在團聚體內部連接較弱的位置更容易發生破裂，破裂后的碎片中顆粒排列比初始團聚體中顆粒排列更為致密．隨著分散的進行，顆粒速度增加，顆粒與流體之間的相對速度減小導致流體拖曳力減小，直至與范德華力達到平衡時不再分散．

圖4 流場強度為820 Pa時團聚體在不同流場中的分散過程Fig．4 Time evolution of dispersion of agglomerates in different flow fields at a flow stress of 820Pa

比較團聚體在不同流場中的分散行為，可以看出，團聚體在拉伸流場中可以獲得較小的碎片，即團聚體更容易分散，混合流場次之，而剪切流場的分散能力最低．這是由于團聚體在剪切流場中受到剪切應力作用后，團聚體分散的同時還存在轉動，消耗了部分流體傳遞給團聚體的能量;此外，團聚體在分散過程中顆粒的排列更為緊密，團聚體的內聚力增加，使其分散變得更為困難．而拉伸流場中團聚體沒有轉動，流體傳遞的能量全部用于團聚體的分散，而且團聚體在分散過程中顆粒排列較為疏松，顆粒間的內聚力較弱，因而團聚體最容易分散．混合流場中的剪切分量亦會使部分能量用于團聚體的轉動，而且分散過程中顆粒的排列也沒有剪切流場中緊密．因此，拉伸流場的分散效率最高，混合流場次之，剪切流場的分散效率最低．

3．2 團聚體分散動力學

為了定量分析流場的分散性能，以分散得到的碎片的加權平均顆粒數〈w〉表征分散度．采用〈w〉隨時間的變化曲線來分析分散混合的動力學過程．定義加權平均顆粒數為

式中，nl為包含l個顆粒的碎片數，l0為團聚體初始時刻的顆粒數．

流場強度為820Pa時，團聚體在剪切、拉伸和剪切-拉伸混合流場中的分散度隨時間的變化如圖5所示．

圖5 流場強度為820Pa時團聚體的分散度隨時間的變化Fig．5 Dispersity of agglomerates changing with time at a flow stress of 820Pa

圖6 不同流場強度時在不同流場下團聚體的分散度隨時間的變化Fig．6 Dispersity of agglomerates changing with time in different flow fields with different flow stresses

由圖5可見:剛開始時團聚體只有變形而沒有分散，經過一定時間后才開始分散;拉伸流場中團聚體最先開始分散，團聚體達到穩定的分散度所用時間最短，分散度也最小;而剪切流場中團聚體從變形到開始分散所用時間最長，達到穩定的分散度所用時間最長，分散度最大．

流場強度為460和1200 Pa時，團聚體在剪切、拉伸和混合流場中的分散度隨時間的變化如圖6所示．

觀察圖6發現，其變化趨勢和圖5一致．比較圖5和6可以看出，對于形式相同的流場，隨著流場強度的增加，由于流場作用于團聚體的能量增加，團聚體從變形到開始分散所用時間縮短，達到穩定的分散度所用時間縮短，分散度減?。?/p>

4 結論

通過對二維情況下均質納米級球形顆粒組成的具有分形結構的團聚體在純剪切、純拉伸和剪切-拉伸混合流場中的分散混合行為進行數值模擬，得出以下結論:

(1)分形結構團聚體在流場中先發生變形然后再分散，分散過程中存在剝蝕和破裂兩種分散方式，分散后的碎片比初始團聚體更加致密．

(2)流場強度相同時，團聚體在拉伸流場中的分散效率最高，碎片平均尺寸最小;在剪切流場中的分散效率最低，碎片平均尺寸最大．

(3)在相同形式的流場中，隨著流場強度增加，團聚體的分散效率提高，碎片平均尺寸減?。?/p>

［1］塔德莫爾Z，高戈斯C G．聚合物加工原理［M］．任冬云，譯．北京:化學工業出版社，1990:239-241．

［2］Feke D L，Manas-Zloczower I．Rupture of inhomogeneous spherical clusters by simple flows［J］．Chemical Engineering Science，1991，46(8):2153-2156．

［3］Horwatt S W，Manas-Zloczower I，Feke D L．Dispersion behavior of heterogeneous agglomerates at supercritical stresses［J］．Polymer Engineering and Science，1992，47(8):1849-1855．

［4］Scurati A，Feke D L，Manas-Zloczower I．Analysis of the kinetics of agglomerate erosion in simple shear flows［J］．Chemical EngineeringScience，2005，60(23):6564-6573．

［5］孫其誠，王光謙．顆粒流動力學及其離散模型評述［J］．力學進展，2008，38(1):87-100．Sun Qi-cheng，Wang Guang-qian．Review on granular flow dynamics and its discrete element method ［J］．Advances in Mechanics，2008，38(1):87-100．

［6］Cundall P A．A computer model for simulating progressive，large-scale movements in blocky rock systems［C］∥Proceedings of the Symposium of the International Society of Rock Mechanics．Rotterdam:Balkama A A，1971:8-12．

［7］徐泳，孫其誠，張凌，等．顆粒離散元法研究進展［J］．力學進展，2003，33(2):251-260．Xu Yong，Sun Qi-cheng，Zhang Ling，et al．Advances in discrete element methods for particulate materials［J］．Advances in Mechanics，2003，33(2):251-260．

［8］Zhu H P，Zhou Z Y，Yang R Y，et al．Discrete particle simulation of particulate systems:theoretical developments［J］．Chemical Engineering Science，2007，62(23):3378-3396．

［9］Zhu H P，Zhou Z Y，Yang R Y，et al．Discrete particle simulation of particulate systems:a review of major applications and findings［J］．Chemical Engineering Science，2008，63(23):5728-5770．

［10］Higashitani K，Iimura K．Two-dimensional simulation of the breakup process of aggregates in shear and elongational flows［J］．Journal of Colloid and Interface Science，1998，204(2):320-327．

［11］Eggersdorfer M L，Kadau D，Herrmann H J，et al．Fragmentation and restructuring of soft-agglomerates under shear［J］．Journal of Colloid and Interface Science，2010，342(2):261-268．

［12］Hosseini M S，Nazockdast H，Dabi B．Numerical simulation of aggregate dispersion in different flow fields using discrete element method ［J］．Applied Polymer Science，2010，115(6):3303-3310．

［13］Cundall P A，Strack O D L．A discrete numerical model for granular assemblies ［J］．Geotechnique，1979，29(1):47-65．

［14］Yang R Y，Zou R P，Yu A B．Computer simulation of the packing of fine particles［J］．Physical Review E，2000，62(3):3900-3908．

［15］林師沛，趙洪，劉芳．塑料加工流變學及其應用［M］．北京:國防工業出版社，2008:58-62．

［16］孫其誠，王光謙．顆粒物質力學導論［M］．北京:科學出版社，2009:63-70．