蘇曉輝，張弛，徐志鋒，金輝，王治國(guó)

（1 西安石油大學(xué)機(jī)械工程學(xué)院，陜西西安 710065； 2 西安交通大學(xué)動(dòng)力工程及多相流國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西西安 710049）

引 言

黏彈性流體懸浮顆粒流廣泛存在于自然界和工業(yè)生產(chǎn)制造過程，如生物體液中細(xì)胞運(yùn)動(dòng)[1-2]、聚合物和日化用品加工、油氣藏水力壓裂[3-4]、微流控裝置中顆粒的連續(xù)操控與分離[5-8]。表面活性劑溶液達(dá)到一定濃度時(shí)表現(xiàn)出很強(qiáng)的黏彈性，具有獨(dú)特的流變學(xué)性質(zhì)，在工業(yè)生產(chǎn)和制造領(lǐng)域受到廣泛關(guān)注。研究人員[9-10]發(fā)現(xiàn)多元聚合物-表面活性劑中表面活性劑的用量對(duì)于原油乳狀液的穩(wěn)定性具有重要影響，溶液彈性能夠有效提高驅(qū)油效率。此外，表面活性劑溶液具有湍流減阻、強(qiáng)化傳熱作用，在流體管道輸送和換熱領(lǐng)域也得到了廣泛應(yīng)用[11-12]。由于其強(qiáng)非線性應(yīng)力應(yīng)變率特性和應(yīng)力松弛效應(yīng)，表現(xiàn)出十分復(fù)雜的流動(dòng)現(xiàn)象[11]。當(dāng)黏彈性表面活性劑溶液中有顆粒存在時(shí)，體系的流動(dòng)特性更為復(fù)雜，例如顆粒在其中沉降時(shí)顆粒速度的“過沖”現(xiàn)象[13]，顆粒沉降速度的持續(xù)波動(dòng)[14-16]等。目前對(duì)于具有復(fù)雜流變性質(zhì)的黏彈性表面活性劑溶液中顆粒沉降規(guī)律仍缺乏清晰的認(rèn)識(shí)，顆粒所呈現(xiàn)的復(fù)雜沉降現(xiàn)象背后的機(jī)理尚不明確。

黏彈性流體中顆粒沉降特性取決于流體的流變性質(zhì)以及操作條件。國(guó)內(nèi)外學(xué)者[17-23]針對(duì)黏彈性流體中顆粒運(yùn)動(dòng)特性開展了相關(guān)研究，主要集中于剪切條件下流體中顆粒橫向運(yùn)移以及成鏈特性，探索了顆粒宏觀運(yùn)動(dòng)與流體流變性質(zhì)之間的關(guān)系。顆粒群在流體中運(yùn)動(dòng)主要取決于流體作用于顆粒的曳力，為此學(xué)者們[24-25]對(duì)黏彈性溶液中顆粒曳力進(jìn)行了研究，發(fā)現(xiàn)彈性(用Deborah 數(shù)De=λu/a來表征，其中λ為松弛時(shí)間；u為流體特征速度；a為特征尺寸，取顆粒半徑)對(duì)顆粒曳力具有顯著影響，在常黏度的Boger 流體中，曳力隨著De增加先減小至某一最小值后又增加，而在剪切稀化流體中，曳力隨著De單調(diào)減小，并且減小速度和幅度與剪切稀化的強(qiáng)度相關(guān)?；诖?，提出了黏彈性流體中顆粒曳力理論公式，用于描述兩相流動(dòng)相間動(dòng)量交換[24,26]。此外，對(duì)黏彈性流體中顆粒沉降特性開展了研究[2,27]。有學(xué)者[13,28-30]采用虛擬區(qū)域法研究黏彈性流體中單顆粒和多顆粒的沉降特性，發(fā)現(xiàn)流體彈性和壁面效應(yīng)均會(huì)影響顆粒沉降，他們均采用只考慮黏彈性的Oldroyd-B 本構(gòu)模型。Oldroyd-B 本構(gòu)模型由于可以簡(jiǎn)潔明確地描述流體的黏彈性，在涉及描述表面活性劑溶液體系的研究中得到了廣泛應(yīng)用[10,13]。

然而，實(shí)際工程中的黏彈性表面活性劑溶液，不僅具有時(shí)間松弛效應(yīng)，其應(yīng)力應(yīng)變率響應(yīng)也具有非線性特性，如剪切稀化、拉伸硬化和剪切帶等，其對(duì)顆粒沉降也有顯著影響。學(xué)者們[14-15,31]通過實(shí)驗(yàn)對(duì)一定濃度下的表面活性劑溶于水后形成的蠕蟲狀膠束溶液中顆粒沉降進(jìn)行研究，發(fā)現(xiàn)顆粒在這種黏彈性溶液中沉降速度呈現(xiàn)隨機(jī)波動(dòng)，不能達(dá)到穩(wěn)定的終端速度，并且認(rèn)為這一現(xiàn)象與顆粒尾跡區(qū)的復(fù)雜拉伸應(yīng)變特性相關(guān)。Wang 等[16]也在實(shí)驗(yàn)中發(fā)現(xiàn)了這一效應(yīng)，并且發(fā)現(xiàn)在某些條件下，不僅顆粒速度在沉降過程中隨機(jī)波動(dòng)，顆粒沉降方向也發(fā)生了變化，即顆粒在流體中出現(xiàn)了沉降-反彈的彈跳現(xiàn)象，認(rèn)為這一現(xiàn)象與此類蠕蟲狀膠束流體的流變性質(zhì)密切相關(guān)[3-4,32]，推測(cè)顆粒沉降過程中下游的“負(fù)尾跡”可能是引起顆粒沉降速度波動(dòng)的原因，因此有必要對(duì)其機(jī)理展開研究。

本文采用同時(shí)考慮表面活性劑溶液黏彈性和剪切稀化及拉伸硬化效應(yīng)的FENE-P 模型以及考慮了大分子之間相互作用的用于濃表面活性劑溶液的Giesekus模型，研究顆粒在其中的沉降特性，討論溶液的流變性質(zhì)對(duì)顆粒沉降特性的影響規(guī)律以及顆粒沉降過程中尾跡區(qū)的特征，研究結(jié)果有助于進(jìn)一步認(rèn)識(shí)和解釋黏彈性表面活性劑溶液中顆粒沉降速度的波動(dòng)機(jī)理。

1 數(shù)學(xué)模型

考慮顆粒沉降運(yùn)動(dòng)過程，為了簡(jiǎn)化模擬，選取顆粒為參考坐標(biāo)系，顆粒在沉降過程中受流體曳力以及重力和浮力作用，其運(yùn)動(dòng)由牛頓運(yùn)動(dòng)定律描述。同時(shí)通過非慣性系的體積力將其與描述流體流動(dòng)的偏微分方程耦合。對(duì)于瞬態(tài)、不可壓縮黏彈性表面活性劑溶液的絕熱流動(dòng)，其控制方程包括連續(xù)方程、動(dòng)量守恒方程以及描述附加彈性應(yīng)力的黏彈性本構(gòu)方程。

連續(xù)方程

動(dòng)量方程

其中，ρ為流體密度；u為流體速度；p為流體壓力；μS為溶劑黏度；τP為附加黏彈性應(yīng)力，由黏彈性流體的本構(gòu)方程確定；F為流體在非慣性系中運(yùn)動(dòng)時(shí)的慣性力。對(duì)于小球參考系，慣性力F只有垂直分量Fz

黏彈性本構(gòu)方程

其中，λ為松弛時(shí)間；μP為溶質(zhì)黏度；fr為松弛函數(shù)；fp為黏度系數(shù)。對(duì)于不同的本構(gòu)方程，fr和fp具有不同的形式。

上隨體導(dǎo)數(shù)

Oldroyd-B本構(gòu)方程

Oldroyd-B 模型把大分子溶液看作是牛頓溶劑流體中懸浮的胡克彈簧體系構(gòu)成，由分子動(dòng)力學(xué)理論推導(dǎo)而來，可用來預(yù)測(cè)基本的黏彈性特征，但是其假定黏度為常數(shù)，并且可以無限拉伸，因而在拉伸流動(dòng)中常常得到不符合物理實(shí)際的結(jié)果，是比較簡(jiǎn)單的一種黏彈性模型。

FENE-P本構(gòu)方程

其中，L為伸長(zhǎng)率。

FENE（finitely extensible nonlinear elastic）模型也是基于分子動(dòng)力學(xué)理論，但是將大分子看作是非線性彈簧，并且考慮了大分子的有限拉伸率。基于FENE模型的Peterlin封閉模型同時(shí)考慮了大分子有限拉伸以及流體剪切稀化效應(yīng)。

Giesekus本構(gòu)方程

其中，α為遷移因子。

Giesekus 本構(gòu)模型用于描述半稀和濃溶液，考慮了大分子之間相互作用產(chǎn)生的流動(dòng)阻力，因此方程中增加了二次非線性項(xiàng)。

顆粒運(yùn)動(dòng)由牛頓運(yùn)動(dòng)定律確定

其中，m為顆粒質(zhì)量；z為顆粒位移；Fz為流體對(duì)顆粒的作用力；Fg為顆粒重力。

流體對(duì)顆粒曳力為

其中，I是單位張量；n為球表面的單位法向量；ez為流體流動(dòng)方向單位矢量；?Ωs是球形顆粒的外表面積。

顆粒所受重力為

其中，ρs為顆粒密度；Vs為顆粒體積；a為顆粒半徑。

雷諾數(shù)和維森貝格數(shù)定義為

其中，d為顆粒粒徑；零剪切黏度μ0=μS+μP；特征速度Us取小球的瞬時(shí)沉降速度。

2 數(shù)值模擬方法

模擬流動(dòng)區(qū)域和計(jì)算區(qū)域如圖1所示，半徑為a的球形顆粒位于半徑為R的圓柱體容器的對(duì)稱軸上，顆粒在重力、浮力和曳力作用下從靜止開始沉降，豎直方向上瞬時(shí)沉降速度為Us，初始時(shí)刻t=0 s時(shí)沉降速度Us=0。由于取運(yùn)動(dòng)顆粒自身為參考系，那么對(duì)于計(jì)算域的流體運(yùn)動(dòng)只需考慮流體相對(duì)于顆粒的運(yùn)動(dòng)即可，顆粒相對(duì)于流體具有平均速度Us，因此在遠(yuǎn)離顆粒的計(jì)算域入口處流體相對(duì)于顆粒參考系的速度u=Us，即流體從計(jì)算域下方入口?Ωin以顆粒瞬時(shí)沉降速度Us均勻垂直滿流流入，從上方出口?Ωout流出，模擬中通過改變通道阻塞比a/R表征壁面對(duì)顆粒沉降的影響，考慮到流場(chǎng)的對(duì)稱性，只需求解如圖1所示的二維軸對(duì)稱流動(dòng)，然后通過周向旋轉(zhuǎn)即可得到三維的顆粒在流體中沉降的流場(chǎng)。

圖1 數(shù)值計(jì)算選取的流動(dòng)區(qū)域和計(jì)算區(qū)域Fig.1 The flow domain and computational domain for numerical simulation

采用有限元軟件COMSOL 的弱形式對(duì)控制方程進(jìn)行求解，采用混合Galerkin 方法對(duì)偏微分方程進(jìn)行離散，對(duì)速度分量采用二次形函數(shù)，對(duì)壓力和黏彈性應(yīng)力采用一次形函數(shù)。通道入口和壁面采用恒定速度u=u(U,0,0)和零壓力梯度?p=0，出口采用壓力出口，顆粒表面流體速度u=u(0,0,0)，壓力梯度?p=0。對(duì)于黏彈性應(yīng)力，入口處為均勻速度，出口處充分發(fā)展，黏彈性應(yīng)力分量為零。容器壁和顆粒表面處溶液中表面活性劑分子只對(duì)壁面產(chǎn)生切向應(yīng)力，其法向應(yīng)力滿足n(n·τP)=0，n是邊界的單位法矢量。

3 流體流變參數(shù)擬合及顆粒沉降模擬

采用控制應(yīng)力流變儀(HAKKE-RS300)對(duì)黏彈性表面活性劑溶液[（2.8%(質(zhì)量)OTAC 和0.56%(質(zhì)量) NaSal）]的線性黏彈性和穩(wěn)態(tài)流變進(jìn)行實(shí)驗(yàn)研究，結(jié)果如圖2 所示，由小幅振蕩剪切實(shí)驗(yàn)結(jié)果，可以獲得溶液流變參數(shù)松弛時(shí)間λ，通過穩(wěn)態(tài)剪切實(shí)驗(yàn)結(jié)果可以獲得表征剪切稀化效應(yīng)的參量如Giesekus 模型中遷移因子α和FENE 模型中的拉伸率L，實(shí)驗(yàn)獲得溶液流變參數(shù)的方法詳見文獻(xiàn)[16]。通過實(shí)驗(yàn)發(fā)現(xiàn)流體具有剪切稀化效應(yīng)，根據(jù)線性黏彈性，松弛時(shí)間λ=0.27 s，μ0=1.944 Pa·s，G0=7.285 Pa，采用Giesekus 本構(gòu)模型進(jìn)行參數(shù)擬合，擬合結(jié)果如圖2（a）所示，采用此流變參數(shù)對(duì)顆粒沉降進(jìn)行模擬，顆粒半徑a為0.003 m，容器半徑R為0.031 m，顆粒密度ρs=7800 kg·m-3，對(duì)比不同遷移因子時(shí)的結(jié)果，模擬所得顆粒沉降速度隨時(shí)間的分布如圖2（b）所示，可以看出，表面活性劑溶液中顆粒沉降速度隨時(shí)間的變化明顯不同于牛頓流體，α較小時(shí)（剪切稀化作用較弱），顆粒沉降初期速度出現(xiàn)了明顯的波動(dòng)性，隨著α增加，剪切稀化作用增強(qiáng)，顆粒沉降速度不穩(wěn)定性減弱。

圖2 黏彈性表面活性劑溶液（2.8%OTAC+0.56%NaSal）流變性質(zhì)實(shí)驗(yàn)數(shù)據(jù)及Giesekus本構(gòu)模型擬合結(jié)果(a)；顆粒沉降速度隨時(shí)間的分布(b)Fig.2 The rheological properties from experiments for viscoelastic surfactant solution (2.8%OTAC+0.56%NaSal)and the data fitting using Giesekus constitutive equation(a);The variation of settling velocity of a particle with time(b)

4 模型驗(yàn)證

首先進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證，考慮到球表面附近速度和應(yīng)力梯度較大，在球表面處采用不同密度的網(wǎng)格M1、M2和M3，對(duì)應(yīng)球表面網(wǎng)格數(shù)量分別為106、160 和320，總 網(wǎng) 格 數(shù) 量 分 別 為71818、106238 和112613，圖3 所示為采用不同網(wǎng)格數(shù)量時(shí)顆粒瞬時(shí)沉降速度的計(jì)算結(jié)果，可以看出，所采用的網(wǎng)格均可以很好地捕捉顆粒沉降的物理特征，通過對(duì)比顆粒沉降終速發(fā)現(xiàn)，采用M2網(wǎng)格相對(duì)于M1網(wǎng)格的顆粒沉降終速相對(duì)變化量為0.03%，而M3網(wǎng)格相對(duì)于M2網(wǎng)格的顆粒沉降終速相對(duì)變化量小于0.01%，據(jù)此采用M2網(wǎng)格作為計(jì)算網(wǎng)格，認(rèn)為進(jìn)一步加密網(wǎng)格對(duì)于模擬結(jié)果的影響可以忽略，對(duì)于不同a/R，球表面附近網(wǎng)格數(shù)量始終取160，而對(duì)遠(yuǎn)離顆粒的區(qū)域采用相對(duì)較粗的網(wǎng)格，既能保證計(jì)算結(jié)果的精度，又能顯著減小計(jì)算量。

圖3 網(wǎng)格無關(guān)性驗(yàn)證Fig.3 Grid independence verification

通過本研究模擬結(jié)果與前人研究對(duì)比，驗(yàn)證本文采用的模擬方法的準(zhǔn)確性。對(duì)比結(jié)果如圖4 所示，采用的本構(gòu)方程為Giesekus模型，為了與文獻(xiàn)已有的結(jié)果（采用Oldroyd-B 模型）對(duì)比，假定本構(gòu)方程中遷移因子α為0，此時(shí)，Giesekus 模型退化為Oldroyd-B 模型，可以看出，采用三種不同通道阻塞比a/R為0.1、0.2和0.5時(shí)，模擬獲得的顆粒沉降速度隨時(shí)間的變化與文獻(xiàn)中結(jié)果吻合較好，說明了本研究采用的模擬方法是可靠的。

由圖4可以看出，阻塞比a/R顯著影響顆粒沉降速度，當(dāng)a/R為0.1 時(shí)，顆粒沉降速度在t=0.1 s 處出現(xiàn)一個(gè)明顯的過沖，即沉降速度很快增加到一個(gè)最大值，而后又逐漸減小，在0.4 s 后達(dá)到穩(wěn)定的終端速度0.16 m·s-1。當(dāng)a/R增加到0.2時(shí)，在t=0.1 s處仍然出現(xiàn)了速度過沖，其速度峰值小于a/R=0.1 時(shí)的值，而后速度減小，減小到某個(gè)值后，仍然沒有趨于穩(wěn)定，而是隨著時(shí)間又逐漸增加，在0.6 s 后基本達(dá)到穩(wěn)定的終端速度0.124 m·s-1，其小于a/R=0.1 時(shí)的終端速度值。當(dāng)a/R進(jìn)一步增加到0.5 時(shí)，壁面影響最強(qiáng)，此時(shí)在0.04 s處就出現(xiàn)速度最大值0.10 m·s-1，其最大值相對(duì)于a/R=0.1時(shí)的最大值減小約一半，沉降速度減小到0.022 m·s-1，之后沉降速度又有明顯增加，再次達(dá)到一個(gè)峰值，后又有輕微減小和增加，如此，經(jīng)過兩次波動(dòng)在0.5 s左右才達(dá)到穩(wěn)定的終端速度0.038 m·s-1，其終速值也遠(yuǎn)小于a/R為0.1 和0.2 時(shí)的終速值。由此說明隨著阻塞比a/R增大，即壁面效應(yīng)的增強(qiáng)，顆粒沉降終速減小，而且顆粒沉降速度不穩(wěn)定性增強(qiáng)，這是因?yàn)楸诿嫘?yīng)較強(qiáng)時(shí)，顆粒沉降過程中應(yīng)變率和應(yīng)力較大，其彈性應(yīng)力更大，彈性效應(yīng)表現(xiàn)更為顯著，顆粒沉降速度經(jīng)過波動(dòng)之后，流體中積聚的彈性能量充分耗散后，才達(dá)到穩(wěn)定的終端速度。

圖4 遷移因子α=0時(shí)Giesekus流體中顆粒沉降速度模擬結(jié)果與文獻(xiàn)[13,33-34]對(duì)比Fig.4 Comparison of settling velocity of a particle in Giesekus fluids(mobility factor α=0)between present simulation results and literature data[13,33-34]

5 結(jié)果和討論

5.1 Giesekus 和FENE-P 本構(gòu)模型的流變性質(zhì)模擬

本文采用的Giesekus 模型適用于半稀和濃溶液，F(xiàn)ENE-P 模型適用于稀溶液，兩種模型都考慮了剪切稀化效應(yīng)，為了明晰其流變性質(zhì)對(duì)流體流動(dòng)和顆粒沉降特性的影響，首先對(duì)Giesekus 和FENE-P本構(gòu)模型的流變性質(zhì)進(jìn)行計(jì)算，考慮到所采用本構(gòu)模型的非線性，采用數(shù)值方法求解在平面均勻剪切以及單軸拉伸流動(dòng)下的流體應(yīng)力分布。首先研究了松弛時(shí)間λ對(duì)流體剪切和拉伸流變性質(zhì)的影響，如圖5所示，可以看出，松弛時(shí)間影響流體黏度曲線關(guān)于應(yīng)變率的分布，當(dāng)松弛時(shí)間較大時(shí)，出現(xiàn)剪切稀化和拉伸硬化對(duì)應(yīng)的應(yīng)變率較小，當(dāng)松弛時(shí)間減小時(shí)，出現(xiàn)剪切稀化和拉伸硬化對(duì)應(yīng)的應(yīng)變率增大，說明了彈性對(duì)于應(yīng)力-應(yīng)變率的非線性性質(zhì)也有影響。

圖5 松弛時(shí)間λ對(duì)流變性質(zhì)的影響Fig.5 The effects of the relaxation time λ on rheological properties

此外研究了兩種本構(gòu)模型中參數(shù)對(duì)流變性質(zhì)的影響，如圖6、圖7 所示，此處采用無量綱化坐標(biāo)，對(duì)剪切流動(dòng)，橫坐標(biāo)Wi=λγ˙，γ˙為剪切速率，縱坐標(biāo)τXY為無量綱總剪切應(yīng)力，對(duì)拉伸流動(dòng)，橫坐標(biāo)Wie=λε˙，ε˙為拉伸速率，縱坐標(biāo)τXX為無量綱總拉伸應(yīng)力，其中λ=1.0 s，μS/μ0=0.1?？梢钥闯?，對(duì)于兩種本構(gòu)模型，在剪切速率0.1 s-1以上，都呈現(xiàn)出蠕蟲狀膠束溶液的典型的剪切稀化效應(yīng)，并且同樣剪切速率，Giesekus 剪切稀化效應(yīng)更強(qiáng)，遷移因子α在0.8 以上時(shí)，Giesekus 流體出現(xiàn)剪切帶。此外，對(duì)于拉伸流動(dòng)，兩種本構(gòu)模型均呈現(xiàn)拉伸硬化的特征，如圖6(b)和圖7(b)所示。

圖6 Giesekus流體在平面均勻剪切(a)和單軸拉伸流動(dòng)(b)時(shí)的流變性質(zhì)（Wi=λγ˙，Wie=λε˙）Fig.6 The dimensionless stress-strain relationships for Giesekus fluids under the planar shear (a)and the uniaxial tension(b)condition (Wi=λγ˙，Wie=λε˙)

圖7 FENE-P流體在平面均勻剪切(a)和單軸拉伸流動(dòng)(b)時(shí)的流變性質(zhì)（Wi=λγ˙，Wie=λε˙）Fig.7 The dimensionless stress-strain relationships for FENE-P fluids under the planar shear(a)and the uniaxial tension(b)condition(Wi=λγ˙，Wie=λε˙)

5.2 顆粒沉降特性

為了闡明流體流變性質(zhì)對(duì)沉降的影響，先不考慮壁面影響，即a/R趨近于0，為此取圓柱通道的半徑R=200a，計(jì)算表明，當(dāng)R/a從150增加到200時(shí)，模擬所得顆粒曳力系數(shù)只變化了0.1%，據(jù)此認(rèn)為繼續(xù)增加R/a對(duì)計(jì)算結(jié)果影響可以忽略不計(jì)，因此取R/a=200，此時(shí)可以看作顆粒在無界流場(chǎng)中沉降，無壁面效應(yīng)影響。研究顆粒浸沒在表面活性劑溶液中運(yùn)動(dòng)特性時(shí)，認(rèn)為顆粒為完全潤(rùn)濕[13-16]。此外，需要指出對(duì)于實(shí)際的顆粒懸浮體系，顆粒群效應(yīng)和顆粒形狀的不規(guī)則均會(huì)對(duì)顆粒沉降特性產(chǎn)生影響，為了簡(jiǎn)化模型，本研究考慮簡(jiǎn)單球形顆粒的沉降來解釋顆粒沉降的非穩(wěn)態(tài)效應(yīng)。

5.2.1 流變參數(shù)對(duì)顆粒沉降速度的影響 在剔除了容器壁面效應(yīng)的影響后，顆粒沉降速度只受表面活性劑溶液流變性質(zhì)的影響，因此研究不同De、α以及L對(duì)顆粒沉降的影響。

先不考慮剪切稀化效應(yīng)，即流體本構(gòu)模型退化為Oldroyd-B 方程，此時(shí)顆粒沉降速度隨時(shí)間的變化規(guī)律如圖8 所示，取橫軸為無量綱時(shí)間因子t/λ，縱軸為無量綱速度u/Ut（Ut為斯托克斯沉降速率），當(dāng)松弛時(shí)間λ為0.01 s 時(shí)，顆粒整體沉降速度較小，速度的增加較為緩慢，經(jīng)過較長(zhǎng)的時(shí)間才能達(dá)到一個(gè)穩(wěn)定的終端速度，隨著λ增加，顆粒初期沉降速度逐漸增大，這是因?yàn)楫?dāng)顆粒沉降過程中對(duì)流體施加應(yīng)變時(shí)，由于應(yīng)力松弛效應(yīng)黏彈性應(yīng)力的增加有一個(gè)滯后，此時(shí)流體應(yīng)力對(duì)顆粒作用力較小，因此沉降速度加速度較大，隨著時(shí)間增加，應(yīng)力增加，對(duì)顆粒阻力作用逐漸增加，顆粒加速度逐漸減小，初始階段的應(yīng)力松弛效應(yīng)，導(dǎo)致了流體中沉降時(shí)顆粒速度的過沖現(xiàn)象，因此彈性效應(yīng)是造成顆粒沉降不穩(wěn)定性的原因。λ增加至0.05 s時(shí)，顆粒沉降速度出現(xiàn)了輕微的過沖現(xiàn)象，當(dāng)λ繼續(xù)增加時(shí)，顆粒沉降速度過沖幅度顯著增加，并且波動(dòng)周期數(shù)也增加，這在壁面約束時(shí)的顆粒沉降速度分布（圖4）上更加明顯。由無量綱時(shí)間看出，隨著λ增加，過沖量達(dá)到最大值的無量綱時(shí)間越小，達(dá)到終速的無量綱時(shí)間越小，這是由于隨著λ增加，流體彈性應(yīng)力松弛效應(yīng)增強(qiáng)，顆粒沉降過程中，流體發(fā)生變形后相應(yīng)的變形應(yīng)力增加緩慢，作用于顆粒的曳力增加緩慢，導(dǎo)致顆粒初始加速度和速度較大，因此產(chǎn)生較為明顯的速度過沖。這說明顆粒在沉降初期的不穩(wěn)定性主要是由流體自身的彈性效應(yīng)引起，彈性效應(yīng)越強(qiáng)，顆粒沉降速度隨時(shí)間不穩(wěn)定性越強(qiáng)。

圖8 無剪切稀化時(shí)松弛時(shí)間λ對(duì)顆粒沉降瞬時(shí)速度的影響Fig.8 The effects of relaxation time λ on transient settling velocity of a particle without shear thinning

除了彈性效應(yīng)，對(duì)剪切稀化對(duì)顆粒沉降的影響進(jìn)行了研究，隨著Giesekus 本構(gòu)方程遷移因子α的增加，剪切稀化效應(yīng)增強(qiáng)，由圖9 可以看出，隨著α增加，剪切稀化增強(qiáng)會(huì)削弱顆粒沉降速度的波動(dòng)，當(dāng)α大于0.3 時(shí)，顆粒沉降速度不再出現(xiàn)過沖現(xiàn)象，而是緩慢增加最后達(dá)到穩(wěn)定的終端速度，這是由于剪切稀化效應(yīng)導(dǎo)致隨著顆粒沉降速度增加，顆粒周圍剪切速率增大，表觀黏度逐漸下降，彈性能量釋放后的應(yīng)力也會(huì)減小，因此顆粒會(huì)有一個(gè)逐漸加速的過程，剪切稀化效應(yīng)越強(qiáng)，顆粒沉降終速越大。

圖9 Giesekus流體遷移因子α對(duì)顆粒沉降瞬時(shí)速度的影響Fig.9 The effects of mobility factor α on transient settling velocity of a particle in Giesekus fluids

對(duì)FENE-P 流體，隨著伸長(zhǎng)率L的增加，剪切稀化效應(yīng)減弱。由圖10 可以看出，拉伸率越小，顆粒沉降速度波動(dòng)幅度越小，當(dāng)L=2時(shí)，顆粒速度不出現(xiàn)過沖現(xiàn)象，因此彈性應(yīng)力松弛效應(yīng)增強(qiáng)顆粒沉降速度的波動(dòng)，而剪切稀化削弱顆粒沉降速度的波動(dòng)。

圖10 FENE-P流體拉伸率L對(duì)顆粒沉降瞬時(shí)速度的影響Fig.10 The effects of extensibility L on transient settling velocity of a particle in FENE-P fluids

在實(shí)驗(yàn)中觀察到顆粒速度的持續(xù)擬周期性波動(dòng)，可以推斷顆粒速度波動(dòng)是由流體的彈性應(yīng)力松弛效應(yīng)觸發(fā)，彈性越強(qiáng)，顆粒沉降波動(dòng)效應(yīng)越強(qiáng)，而剪切稀化效應(yīng)會(huì)削弱顆粒的波動(dòng)，但是采用Giesekus 和FENE-P 本構(gòu)方程模擬的結(jié)果發(fā)現(xiàn)顆粒在開始沉降后，速度波動(dòng)很快消失，達(dá)到穩(wěn)定的終端速度。黏彈性流動(dòng)的高維森貝格數(shù)會(huì)引起求解不穩(wěn)定，容易發(fā)散，本文在模擬中發(fā)現(xiàn)De超過10 后計(jì)算變得難以收斂，使得模擬更強(qiáng)彈性效應(yīng)的沉降變得困難，可以推斷，隨著彈性進(jìn)一步增加，顆粒沉降速度波動(dòng)周期數(shù)和幅度均增加。另一方面，表面活性劑分子形成的蠕蟲狀膠束溶液由于其流變性質(zhì)比共價(jià)骨架連接的聚合物分子更為復(fù)雜，其是一種活性聚合物，膠束之間不斷地?cái)嗔押椭貥?gòu)，其斷裂-重構(gòu)效應(yīng)會(huì)對(duì)顆粒運(yùn)動(dòng)產(chǎn)生不可忽略的影響，需要從流體中表面活性劑溶液形成的微觀膠束結(jié)構(gòu)的動(dòng)態(tài)演變規(guī)律去揭示彈性和膠束斷裂-重構(gòu)的相互作用機(jī)理。

5.2.2 顆粒沉降負(fù)尾跡特性 表面活性劑溶液中顆粒沉降速度的波動(dòng)可能與顆粒負(fù)尾跡有關(guān)，因此對(duì)顆粒負(fù)尾跡特性進(jìn)行了研究，為了獲得顆粒負(fù)尾跡的流場(chǎng)特征，將模擬結(jié)果轉(zhuǎn)換到慣性參考系，即靜止流體中小球的自然下落，結(jié)果如圖11 所示，可以看出，顆粒在沉降過程中其后緣滯止點(diǎn)（相對(duì)于顆粒速度為零的點(diǎn)）附近的流體由于黏性作用隨著顆粒向下運(yùn)動(dòng)，前緣滯止點(diǎn)附近的流體被顆粒推動(dòng)向下運(yùn)動(dòng)，但是在顆粒兩側(cè)表面附近以及向下游延伸的區(qū)域里，流體速度反而向上，這是由于此區(qū)域的流體受剪切作用較強(qiáng)，由于流體的剪切稀化效應(yīng)，顆粒周圍向下的流體對(duì)周圍流體的黏滯作用力減小，因此速度最大的區(qū)域?qū)ν鈬挠绊憛^(qū)域也減小，強(qiáng)剪切帶把此區(qū)域分成兩個(gè)部分，就像顆粒帶動(dòng)尾跡區(qū)將流體劈開一樣。此外，在尾跡區(qū)向下游延伸處，出現(xiàn)了一個(gè)速度方向與顆粒沉降方向相反的區(qū)域，如圖11(a)中“負(fù)尾跡”所示，這部分流體不斷地卷入尾跡區(qū)中的流體向上運(yùn)動(dòng)，經(jīng)過一段距離后又被主流攜帶進(jìn)入負(fù)尾跡旁邊的主流區(qū)，轉(zhuǎn)而速度向下，在顆粒下游處形成一個(gè)循環(huán)區(qū)域，形成負(fù)尾跡區(qū)，此處對(duì)附近流體的剪切和拉伸會(huì)造成大量膠束的斷裂和重構(gòu)，使得流體內(nèi)部的應(yīng)力變得不穩(wěn)定，是造成流場(chǎng)不穩(wěn)定及顆粒沉降速度波動(dòng)的根本原因。圖11(b)、(c)為沿對(duì)稱軸的顆粒尾跡區(qū)中心線上流體速度分布，縱坐標(biāo)為速度大小，橫坐標(biāo)為從顆粒后緣滯止點(diǎn)開始的距離，在沉降過程中顆粒表面附近流體由于流體黏滯作用隨著顆粒一起向下運(yùn)動(dòng)，速度為負(fù)，在某些區(qū)域出現(xiàn)速度大于零，流體速度與顆粒沉降速度相反，即為負(fù)尾跡。由圖11(b)可以看出，隨著α增大，剪切稀化和拉伸硬化效應(yīng)增強(qiáng)，負(fù)尾跡區(qū)增大，并且距離顆粒更遠(yuǎn)。對(duì)Giesekus 流體，當(dāng)α為0.1 時(shí)，負(fù)尾跡區(qū)反向速度值最大，但是負(fù)尾跡向下游延伸較短，隨著α逐漸增加至0.7，反向速度值減小，負(fù)尾跡向下游延伸加長(zhǎng)，隨著α繼續(xù)增加，負(fù)尾跡反向速度又逐漸增加。對(duì)FENE-P 流體，隨著拉伸率減小，剪切稀化增強(qiáng)，反向速度峰值和尾跡區(qū)向下延伸長(zhǎng)度也是先增大后減小[圖11(c)]。

圖11 以速度矢量表示的顆粒下游負(fù)尾跡(a)；Giesekus模型遷移因子α對(duì)負(fù)尾跡的影響(b)；FENE-P模型拉伸率L對(duì)負(fù)尾跡的影響(c)Fig.11 The negative wake downstream of the particle presented using velocity vector(a);The effects of mobility factor α on negative wake for Giesekus fluids(b);The effects of extensibility L on negative wake for FENE-P fluids(c)

因此，流體剪切稀化和拉伸硬化共同導(dǎo)致顆粒周圍流場(chǎng)出現(xiàn)負(fù)尾跡，其對(duì)顆粒沉降速度不穩(wěn)定性有著重要作用。

此外，研究了在同樣剪切稀化和拉伸硬化時(shí)彈性對(duì)負(fù)尾跡的影響，如圖12 所示，可以看出，對(duì)于同樣的剪切稀化和拉伸硬化（對(duì)Giesekus 流體遷移因子α取0.1，對(duì)FENE-P 流體伸長(zhǎng)率L取5），當(dāng)松弛時(shí)間增加，即彈性效應(yīng)增強(qiáng)時(shí)，顆粒負(fù)尾跡反向速度峰值和負(fù)尾跡向下游延伸長(zhǎng)度都增加，說明了彈性效應(yīng)越強(qiáng)，負(fù)尾跡越強(qiáng)，因此流體的高彈性會(huì)強(qiáng)化顆粒負(fù)尾跡，其與剪切稀化和拉伸硬化共同導(dǎo)致了顆粒周圍流體流動(dòng)以及顆粒沉降的不穩(wěn)定，這與剪切稀化和拉伸硬化導(dǎo)致的微觀蠕蟲狀膠束的動(dòng)態(tài)斷裂-重構(gòu)有關(guān)，為完全準(zhǔn)確描述蠕蟲狀膠束溶液中顆粒沉降的波動(dòng)現(xiàn)象，需要進(jìn)一步探討適用于微觀膠束動(dòng)態(tài)斷裂-重構(gòu)機(jī)制的本構(gòu)模型。

圖12 黏彈性對(duì)顆粒沉降過程中負(fù)尾跡的影響Fig.12 The effects of viscoelasticity on negative wake during the particle settling

6 結(jié) 論

本文對(duì)黏彈性表面活性劑溶液中球形顆粒沉降特性進(jìn)行了研究，同時(shí)考慮了溶液黏彈性效應(yīng)和非線性應(yīng)力應(yīng)變響應(yīng)特性，發(fā)現(xiàn)顆粒在沉降初期的不穩(wěn)定性主要是由流體自身的彈性效應(yīng)引起，彈性效應(yīng)越強(qiáng)，顆粒沉降速度隨時(shí)間不穩(wěn)定性越強(qiáng)，而剪切稀化效應(yīng)會(huì)削弱顆粒沉降速度的波動(dòng)。

此外發(fā)現(xiàn)顆粒在表面活性劑溶液中沉降時(shí)在流體剪切稀化和拉伸硬化作用下顆粒下游形成的“負(fù)尾跡”，而在純黏彈性流體無這一現(xiàn)象發(fā)生，驗(yàn)證了實(shí)驗(yàn)中對(duì)蠕蟲狀表面活性劑溶液中顆粒沉降時(shí)負(fù)尾跡存在的推測(cè)，并且發(fā)現(xiàn)剪切稀化和拉伸硬化效應(yīng)增強(qiáng)會(huì)強(qiáng)化負(fù)尾跡效應(yīng)。彈性效應(yīng)增加時(shí)，負(fù)尾跡區(qū)反向最大速度和向下游延伸長(zhǎng)度均增加，負(fù)尾跡區(qū)與顆粒下沉方向相反的速度分布可能是引起顆粒沉降速度持續(xù)波動(dòng)的原因。

采用Giesekus 和FENE-P 模型雖然可以描述顆粒在經(jīng)歷初始階段的沉降速度的波動(dòng)性，但是在解釋顆粒連續(xù)周期性波動(dòng)時(shí)仍有一定的局限性，這是因?yàn)楸砻婊钚詣┤芤褐行纬傻娜湎x狀膠束是一種“活性聚合物”，不同于高分子聚合物溶液，這種流體在剪切拉伸時(shí)微觀尺度的膠束處于不斷地?cái)嗔?重構(gòu)的動(dòng)態(tài)過程，而傳統(tǒng)的黏彈性本構(gòu)模型不能描述在拉伸過程的斷裂-重構(gòu)特性，因此下一步研究應(yīng)考慮流體中微觀膠束的斷裂-重構(gòu)動(dòng)態(tài)過程對(duì)于顆粒沉降過程的影響，研究能更好地描述表面活性劑蠕蟲狀膠束溶液流變性質(zhì)的本構(gòu)方程。此外，對(duì)于懸浮顆粒體系，顆粒群效應(yīng)和非球形不規(guī)則顆粒均會(huì)對(duì)顆粒沉降和運(yùn)移產(chǎn)生重要影響，因此后續(xù)應(yīng)對(duì)顆粒群效應(yīng)和顆粒形狀的影響開展研究，獲得黏彈性表面活性劑溶液中顆粒沉降的普適規(guī)律。

符 號(hào) 說 明

a——小球半徑，m

D——圓柱通道直徑，m

d——小球直徑，m

ez——z方向單位矢量

F——慣性力，N

Fg——小球重力，N

Fz——慣性力在方向z的分量，N

fp——黏度系數(shù)

fr——松弛函數(shù)

G0——彈性模量，Pa

I——單位張量

L——伸長(zhǎng)率

n——球表面法向矢量

p——壓力，Pa

R——容器半徑，m

Re——雷諾數(shù)

t——時(shí)間，s

Us——顆粒沉降速率，m·s-1

Ut——顆粒斯托克斯沉降速率，m·s-1

u——流體速度矢量，m·s-1

Vs——球形顆粒體積，m3

Wi——維森貝格數(shù)

α——遷移因子

γ˙——流體應(yīng)變率，s-1

λ——松弛時(shí)間，s

μP——溶質(zhì)黏度，Pa·s

μS——溶劑黏度，Pa·s

μ0——零剪切黏度，Pa·s

ρ——流體密度，kg·m-3

ρs——顆粒密度，kg·m-3

τP——黏彈性應(yīng)力，Pa