活性流體流變行為的布朗動(dòng)力學(xué)模擬研究1)

許曉飛 童松豪 張 達(dá) 董 超 劉鳳霞 魏 煒 劉志軍

(大連理工大學(xué)流體與粉體工程研究設(shè)計(jì)所，遼寧大連 116012)

引言

存在于自然界中(如大腸桿菌[1])或是人工合成的(如Janus 粒子[2])活性粒子(active particles)能夠?qū)⑵渲車h(huán)境中各種形式的能量轉(zhuǎn)化為空間運(yùn)動(dòng)的動(dòng)能，驅(qū)動(dòng)自身運(yùn)動(dòng)并擾動(dòng)周圍流場.包含大量活性粒子的懸浮液被稱為活性流體(active fluids).被動(dòng)粒子(passive particles)則是指那些不能自發(fā)運(yùn)動(dòng)，只能通過布朗運(yùn)動(dòng)和流體對流輸運(yùn)發(fā)生遷移的粒子.由于活性粒子的自驅(qū)動(dòng)性，活性流體屬于一種不遵循能量耗散定理的非平衡體系，具有與被動(dòng)流體(passive fluids)完全不同的物理特性，例如生物湍流[3]、相分離行為[4]、超擴(kuò)散現(xiàn)象[5]和特殊流變行為[6]等.深入研究活性流體流變學(xué)不僅有助于理解生物物理學(xué)的基本過程、更好地控制微生物輸運(yùn)過程中生物膜的形成以及改善微流控設(shè)備中的混合過程，更重要的是可以為當(dāng)前新材料的開發(fā)提供設(shè)計(jì)幫助[7-8].活性粒子一般為亞微米級或微米級[9-10]，當(dāng)前實(shí)驗(yàn)研究主要專注于宏觀流變行為[11-12]，而借助數(shù)值模擬的方法可以通過再現(xiàn)活性粒子的自驅(qū)動(dòng)行為，從介觀尺度分析宏觀流變行為的形成機(jī)理，也能研究外部剪切流與活性粒子自驅(qū)動(dòng)行為之間的耦合作用，因此逐漸引起學(xué)者們的關(guān)注.

活性粒子可以響應(yīng)外部流場的剪切作用，在活性流體中產(chǎn)生“潤滑”作用，使懸浮液的黏度下降，甚至有可能使宏觀黏度降至液相黏度之下[13-16].Sokolov和Aranson[17]首次完成細(xì)菌懸浮液流變學(xué)的測試，結(jié)果發(fā)現(xiàn)濃細(xì)菌懸浮液的黏度甚至低于未加入細(xì)菌的基液的黏度.Gachelin 等[18]利用Y 型微通道實(shí)驗(yàn)研究了25 °C 下大腸桿菌懸浮液的相對黏度與剪切速率之間的關(guān)系，結(jié)果發(fā)現(xiàn)在低剪切速率(＜10 s-1)條件下，大腸桿菌懸浮液的黏度低于液相黏度，且隨著剪切速率的增大，懸浮液的相對黏度不斷增大，即出現(xiàn)剪切增稠現(xiàn)象，在臨界剪切速率20 s-1處存在最大值;當(dāng)剪切速率大于臨界剪切速率時(shí)又出現(xiàn)了剪切稀化現(xiàn)象.

為了揭示活性流體特殊流變行為的成因，數(shù)值模擬的手段被廣泛采用以研究活性粒子介尺度行為對體系宏觀流變行為的影響.Saintillan 和Shelley[19]基于細(xì)長體理論創(chuàng)建了活性流體模型，結(jié)果發(fā)現(xiàn)該體系具有集群運(yùn)動(dòng)的特性.Ezhilan 和Saintillan[20]建立了控制懸浮液中粒子速度、濃度和應(yīng)力場的控制方程，該控制方程適用于細(xì)菌懸浮液等一系列活性流體.Haines 等[21]建立活性流體的二維模型，引入了粒子自驅(qū)動(dòng)項(xiàng)，模擬的結(jié)果與實(shí)驗(yàn)類似，觀察到了活性流體的黏度下降現(xiàn)象.后來他們還將模型推廣到三維[22]，考慮了活性粒子的轉(zhuǎn)向行為以及由于粒子形狀不對稱引起的隨機(jī)扭矩的影響，結(jié)果表明活性粒子的自驅(qū)動(dòng)行為降低了有效黏度.Ryan 等[23]以偶極子的方式創(chuàng)建黏性流體中的稀釋活性粒子體系，不考慮粒子間的相互作用.模擬結(jié)果表明，隨著活性粒子濃度的增加，懸浮體系的黏度開始逐漸減小，達(dá)到一個(gè)臨界點(diǎn)后開始增加，這一變化趨勢與Sokolov 和Aranson[17]的實(shí)驗(yàn)結(jié)果一致;然而，這種直接忽略短程相互作用，僅考慮粒子與流體的相互作用的方式可能導(dǎo)致模擬結(jié)果的偏差.Pagonabarraga和Llopis[24]模擬庫埃特流中偶極蠕動(dòng)粒子懸浮液和無極蠕動(dòng)粒子懸浮液的穩(wěn)態(tài)剪切響應(yīng)，結(jié)果表明前者出現(xiàn)剪切稀化特性，后者出現(xiàn)剪切增稠特性，認(rèn)為活性流體黏度下降是粒子內(nèi)部的活性力與流體阻力競爭的結(jié)果，然而并未直接給出黏度的結(jié)果，并且模擬結(jié)果與流變學(xué)實(shí)驗(yàn)結(jié)果不一致.

過往的模擬研究為了簡化計(jì)算，通常會(huì)忽略一些影響因素，例如粒子間作用力、活性粒子的轉(zhuǎn)向運(yùn)動(dòng)和直行速度的變化等，因此會(huì)出現(xiàn)和實(shí)驗(yàn)結(jié)果不一致的情況.此外，還存在不能發(fā)揮模擬計(jì)算單因素分析優(yōu)勢的問題，僅考慮活性粒子體積分?jǐn)?shù)的影響，而忽略了粒子自驅(qū)動(dòng)的影響.本文基于粗粒化的布朗動(dòng)力學(xué)方法[25]，在布朗運(yùn)動(dòng)方程中添加直行運(yùn)動(dòng)項(xiàng)并耦合粒子轉(zhuǎn)向運(yùn)動(dòng)，建立可以描述活性粒子運(yùn)動(dòng)的活性布朗動(dòng)力學(xué)控制方程;在準(zhǔn)二維模擬區(qū)域中建立活性流體模型，使用反向非平衡法模擬計(jì)算活性流體體系的黏度;考察活性粒子體積分?jǐn)?shù)、直行速度和轉(zhuǎn)向擴(kuò)散系數(shù)對活性流體流變行為的影響規(guī)律，并且分析活性粒子的漲落規(guī)律，探究活性流體特殊流變行為的形成機(jī)理.

1 模型的建立和求解

1.1 布朗運(yùn)動(dòng)方程的擴(kuò)展

對于一個(gè)懸浮在液體中，坐標(biāo)為q，質(zhì)量為m，動(dòng)量為p的剛性球粒子，其運(yùn)動(dòng)可用布朗運(yùn)動(dòng)方程描述

其中，t是時(shí)間;表示粒子受到的摩擦阻力，由該懸浮粒子與周圍溶劑分子摩擦產(chǎn)生(α表示斯托克斯系數(shù)，對于球形粒子，根據(jù)斯托克斯定律有α=6πηa，其中η為流體黏度，a為粒子半徑);ξ(t)=λR(t)表示粒子與流體分子碰撞產(chǎn)生的隨機(jī)力(λ=(2DT)1/2α表示阻尼系數(shù)，其中直行擴(kuò)散系數(shù)，kB是玻爾茲曼常數(shù)，T是體系的絕對溫度;R是高斯隨機(jī)噪聲)，使粒子的運(yùn)動(dòng)軌跡雜亂無章，隨機(jī)力的均值是零，滿足高斯分布，因此也稱為高斯噪聲項(xiàng);除此以外，粒子還會(huì)受到勢力-?U(q)，U是粒子之間的相互作用勢.本文使用Lennard-Jones (LJ)勢函數(shù)描述粒子間相互作用勢

其中，r為粒子間的距離;ε和σ表示勢能參數(shù)，與粒子的種類和性質(zhì)有關(guān).式(2) 中，為排斥項(xiàng)，為吸引項(xiàng).粒子間距離越近，排斥力力越大，通過勢力項(xiàng)-?U(q) 即可反映粒子間碰撞的過程.

活性粒子為亞微米或微米尺寸，其運(yùn)動(dòng)屬于低雷諾數(shù)運(yùn)動(dòng)，因此活性粒子受到的慣性力可以忽略不計(jì).將直行擴(kuò)散系數(shù)DT代入布朗運(yùn)動(dòng)方程中整理得

與被動(dòng)粒子相比，活性粒子在流體中運(yùn)動(dòng)時(shí)會(huì)將系統(tǒng)中其他形式能量轉(zhuǎn)化為自身運(yùn)動(dòng)的動(dòng)能.因此，除了流體分子對粒子的碰撞和摩擦作用以外，活性粒子還在自驅(qū)動(dòng)力的作用下產(chǎn)生一個(gè)直行速度U0n(其中，n表示活性粒子自驅(qū)動(dòng)的方向向量).擴(kuò)展后的布朗運(yùn)動(dòng)方程為

除直行運(yùn)動(dòng)以外，活性粒子的運(yùn)動(dòng)還包括轉(zhuǎn)向運(yùn)動(dòng).描述活性粒子轉(zhuǎn)向運(yùn)動(dòng)的方程如下

其中，ψ是活性粒子自驅(qū)動(dòng)方向角，dψ是矢量，角速度ω=dψ/dt也是矢量;是活性粒子特征時(shí)間尺度τR的倒數(shù)，表示粒子的轉(zhuǎn)向擴(kuò)散系數(shù).

活性粒子自驅(qū)動(dòng)方向是隨機(jī)的，活性粒子的實(shí)際運(yùn)動(dòng)是直行運(yùn)動(dòng)和轉(zhuǎn)向運(yùn)動(dòng)的組合.因此，描述活性粒子運(yùn)動(dòng)的方程式為

根據(jù)式(6)建立的數(shù)學(xué)模型，單個(gè)活性粒子的運(yùn)動(dòng)軌跡如圖1 所示，粒子從環(huán)境中獲得能量后所呈現(xiàn)的運(yùn)動(dòng)形式用直行速度U0和轉(zhuǎn)向擴(kuò)散系數(shù)DR來描述.

圖1 活性粒子運(yùn)動(dòng)示意圖Fig.1 Schematic diagram of active particle locomotion

1.2 活性流體黏度的模擬計(jì)算

分子動(dòng)力學(xué)模擬計(jì)算黏度的方法主要包括平衡分子動(dòng)力學(xué)法(equilibrium molecular dynamics，EMD)[26]和非平衡分子動(dòng)力學(xué)法(non-equilibrium molecular dynamics，NEMD)[27].其中，EMD 法主要使用Green-Kubo 關(guān)系[28]確定流體黏度

其中，V是系統(tǒng)的體積，Pαβ為應(yīng)力張量的分量.然而，對于活性流體這樣的非平衡體系，使用NEMD 法計(jì)算會(huì)得到更加準(zhǔn)確的結(jié)果.NEMD 法的思路是對流體施加一個(gè)剪切作用，計(jì)算流體在剪切作用下的響應(yīng)，借助計(jì)算動(dòng)量通量和速度梯度的比值確定流體的黏度[29-30].Müller-Plathe[31]對NEMD 法做出改進(jìn)，提出以非物理(unphysical)的方法在系統(tǒng)中施加動(dòng)量通量.Müller-Plathe 的思路是將模擬盒子沿指定方向劃分為M層，第1 層在指定方向的下邊緣，第M層在指定方向的上邊緣.每N步交換模擬盒子中不同層的兩個(gè)粒子的動(dòng)量，由于體系中每個(gè)粒子的質(zhì)量相同，相當(dāng)于二者的動(dòng)能也發(fā)生了交換.隨著模擬過程的進(jìn)行，可以在體系中產(chǎn)生一個(gè)剪切速率曲線來計(jì)算流體的黏度.該方法也被稱為反向非平衡分子動(dòng)力學(xué)(reverse non-equilibrium molecular dynamics，rNEMD)法，這是因?yàn)镹EMD 法是通過對體系施加剪切作用，然后對應(yīng)力張量的非對角元(正比于動(dòng)量通量)進(jìn)行統(tǒng)計(jì)計(jì)算，而rNEMD 法是先施加動(dòng)量通量，剪切速率曲線是流體響應(yīng)的結(jié)果.

在布朗動(dòng)力學(xué)的模擬過程中，通過對流體速度場、特征剪切應(yīng)力和特征剪切速率等物理量的統(tǒng)計(jì)計(jì)算，就可以獲得流體的黏度值.為驗(yàn)證黏度計(jì)算方法的可靠性，采用簡單LJ 流體進(jìn)行測試，其理論黏度η*=1.0.圖2 給出采用EMD 法、NEMD 法和rNEMD 法模擬計(jì)算LJ 流體黏度的結(jié)果.從圖2 可以看出，EMD 法計(jì)算的黏度隨時(shí)間波動(dòng)較大，計(jì)算得到的黏度平均值為0.946;NEMD 法計(jì)算得到的黏度平均值為1.18;而使用反向非平衡的rNEMD 法獲得的被動(dòng)流體黏度平均值為0.997，更接近該流體的理論黏度值，且相對于其他兩種方法計(jì)算結(jié)果更加穩(wěn)定，可以快速準(zhǔn)確地確定體系的黏度.因此，本文采用rNEMD 法計(jì)算活性流體的黏度.

圖2 3 種方法計(jì)算被動(dòng)LJ 流體的黏度隨時(shí)間的變化關(guān)系Fig.2 Variation of viscosity of LJ passive fluid with time simulated by three schemes

1.3 模型的求解計(jì)算

本研究使用LAMMPS 軟件求解描述活性粒子運(yùn)動(dòng)的擴(kuò)展布朗運(yùn)動(dòng)方程，采用Velocity-Verlet 算法，粒子間作用勢能選擇LJ 勢函數(shù)，運(yùn)算時(shí)的物理量使用LJ 單位，邊界處采用周期性邊界條件.本文為了考察活性粒子體積分?jǐn)?shù)、直行速度和轉(zhuǎn)向擴(kuò)散系數(shù)對活性流體流變行為的影響規(guī)律，確定模擬區(qū)域?yàn)?.5 mm × 0.5 mm × 0.1 mm 的準(zhǔn)二維模擬盒，活性粒子被設(shè)置為剛性球體，活性流體中粒子的體積分?jǐn)?shù)φ分別為0.05%，0.1%，0.3%和0.5%，活性粒子的直行速度U0的大小分別為0 μm/s (無直行運(yùn)動(dòng))、20 μm/s，30 μm/s 和40 μm/s，活性粒子的轉(zhuǎn)向行為由轉(zhuǎn)向擴(kuò)散系數(shù)DR反映，大小分別為0 s-1(無轉(zhuǎn)向運(yùn)動(dòng))、1.0 × 10-4s-1、5.0 × 10-4s-1和1.0 × 10-3s-1.

通過求解描述活性粒子運(yùn)動(dòng)的擴(kuò)展布朗運(yùn)動(dòng)方程，可以得到體系中每個(gè)活性粒子的運(yùn)動(dòng)軌跡信息.活性流體的非平衡特性和活性布朗動(dòng)力學(xué)中的隨機(jī)機(jī)制使得結(jié)果的準(zhǔn)確性受到影響，需要大量的重復(fù)統(tǒng)計(jì)以獲取更為接近實(shí)際情況的結(jié)果.因此，本文每一種工況的模擬計(jì)算，先運(yùn)行5.0 × 104τ(τ是體系特征時(shí)間)的時(shí)間進(jìn)行弛豫，然后運(yùn)行2.0 × 105τ的時(shí)間以獲取所需的宏觀物理信息，計(jì)算中統(tǒng)計(jì)時(shí)間步長Δt=0.001τ.使用大量獨(dú)立軌跡的時(shí)間平均來獲得系統(tǒng)平均，這對于體系宏觀物理量的計(jì)算是必要的.

1.4 模型的可靠性驗(yàn)證

被動(dòng)流體是平衡體系，高斯噪聲的均值為零，體系內(nèi)粒子分布是均勻的，只有在外部輸入能量的情況下才會(huì)出現(xiàn)相變.而活性流體屬于非平衡體系，無需外界輸入能量，通過粒子間的相互作用就可發(fā)生“相變”，出現(xiàn)分離的“液相”和“氣相”，具體表現(xiàn)為粒子聚集，并出現(xiàn)分布不均勻的情況.圖3 給出活性流體和被動(dòng)流體中粒子分布的對比.從圖3 可以看出，被動(dòng)流體作為平衡體系，粒子的運(yùn)動(dòng)遵從布朗運(yùn)動(dòng)方程，長時(shí)間運(yùn)動(dòng)后仍呈現(xiàn)均勻分布狀態(tài)(見圖3(b));而活性流體是非平衡體系，會(huì)自發(fā)地出現(xiàn)粒子聚集的“液相”區(qū)和粒子稀疏的“氣相”區(qū)(見圖3(c))，再現(xiàn)非平衡體系的這一特征證明了本文提出的擴(kuò)展布朗運(yùn)動(dòng)方程的正確性.

圖3 被動(dòng)流體平衡體系和活性流體非平衡體系粒子分布的對比Fig.3 Comparison of particle distribution between equilibrium system of passive fluids and non-equilibrium system of active fluids

因易于配制、穩(wěn)定性高和粒子自驅(qū)動(dòng)機(jī)制簡單等特點(diǎn)，大腸桿菌懸浮液通常被作為典型的活性流體模型來研究活性流體所表現(xiàn)出的一些特殊的物理現(xiàn)象.基于此，大腸桿菌懸浮液的流變學(xué)實(shí)驗(yàn)測試結(jié)果被用來驗(yàn)證本文提出的數(shù)學(xué)模型和模擬計(jì)算方法的可靠性.設(shè)不加入活性粒子的液相黏度為η0，加入活性粒子后體系的黏度為η，定義體系的相對黏度為ηr=η/η0.圖4 給出使用rNEMD 法計(jì)算的活性流體相對黏度結(jié)果與López 等[13]的實(shí)驗(yàn)結(jié)果的對比.自驅(qū)動(dòng)活性粒子的直行速度U0和轉(zhuǎn)向擴(kuò)散率D均使用López 等[13]實(shí)驗(yàn)中測得的數(shù)據(jù)，其中U0=20 μm/s，D=13 μm2/s.從圖4 可以看出，在低剪切速率區(qū)域，基于本文確定的活性流體模型及rNEMD 法所得的結(jié)果與實(shí)驗(yàn)結(jié)果一致性良好.與實(shí)驗(yàn)手段不同的是，布朗動(dòng)力學(xué)模擬的方法可以更容易地控制變量，分析多因素對活性流體流變行為的影響規(guī)律.

圖4 相對黏度隨剪切速率變化的模擬結(jié)果和實(shí)驗(yàn)數(shù)據(jù)對比Fig.4 Comparison of simulated results and experimental data of relative viscosity of active fluids

2 模擬結(jié)果與分析

2.1 活性流體流變行為的影響因素分析

圖5 給出不同體積分?jǐn)?shù)活性流體的相對黏度ηr隨剪切速率的變化.活性粒子的自驅(qū)動(dòng)速度U0=40 μm/s，轉(zhuǎn)向擴(kuò)散系數(shù)DR=0.000 1 s-1.從圖5 可以看出，對于體積分?jǐn)?shù)為0.05%的活性流體，相對黏度隨剪切速率增大沒有明顯改變;而體積分?jǐn)?shù)φ≥0.1%的活性流體呈現(xiàn)非牛頓流體特性，流變曲線具有非單調(diào)性，可以劃分為3 個(gè)區(qū)域，具體包括:(1) 黏度下降區(qū)(＜～ 10 s-1)，低剪切速率區(qū)域內(nèi)活性流體的相對黏度ηr＜1，即活性流體的黏度低于未加入活性粒子的流體黏度，并且活性粒子的體積分?jǐn)?shù)越高，黏度下降越明顯;(2) 過渡區(qū)(～ 10 s-1＜＜～ 50 s-1)，中剪切速率區(qū)域內(nèi)活性流體的黏度出現(xiàn)先增大后減小的變化趨勢，并且存在一個(gè)臨界剪切速率c，此時(shí)活性流體的黏度最大，臨界剪切速率c出現(xiàn)在剪切速率=20 s-1附近;(3) 牛頓區(qū)(＞～ 50 s-1)，高剪切速率區(qū)域內(nèi)活性流體的相對黏度ηr≈ 1，可以被看作為牛頓流體.活性流體的非單調(diào)流變學(xué)行為受到活性粒子自驅(qū)動(dòng)行為和外部施加的剪切應(yīng)力耦合作用的影響，當(dāng)剪切速率較小時(shí)，活性粒子自驅(qū)動(dòng)行為的貢獻(xiàn)占主導(dǎo)，ηr＜ 1;當(dāng)剪切速率較大時(shí)，外部施加的剪切應(yīng)力的作用占主導(dǎo)，ηr≈ 1;當(dāng)二者作用的大小相當(dāng)時(shí)，出現(xiàn)黏度變化的過渡區(qū)域，黏度值小幅度增大.

圖5 活性粒子的體積分?jǐn)?shù)對活性流體相對黏度的影響Fig.5 Effect of volume fraction of active fluid on relative viscosity of active particles

圖6 給出具有不同直行速度活性粒子的活性流體相對黏度隨剪切速率的變化.活性粒子的體積分?jǐn)?shù)φ=0.5%，活性粒子的轉(zhuǎn)向擴(kuò)散系數(shù)DR=0.000 1 s-1.從圖6 可以看出，對于活性粒子直行速度U0=0 μm/s 的活性流體，相對黏度隨剪切速率增大沒有明顯改變;活性粒子的直行速度U0的大小對低剪切區(qū)黏度的下降有顯著貢獻(xiàn)，U0越大則在同一剪切速率下黏度下降的幅度越大;而且活性流體在過渡區(qū)形成的先增稠后稀化現(xiàn)象也更顯著，證明活性粒子的直行運(yùn)動(dòng)對活性流體流變行為有重要影響.

圖6 活性粒子直行速度對活性流體黏度的影響Fig.6 Effect of forward locomotion velocity of active particles on relative viscosity of active fluids

除了直行速度外，控制其自驅(qū)動(dòng)方向改變頻率的轉(zhuǎn)向擴(kuò)散系數(shù)DR同樣是描述活性粒子自驅(qū)動(dòng)行為的重要參數(shù).當(dāng)DR很大時(shí)，活性粒子會(huì)在原地不斷轉(zhuǎn)向，導(dǎo)致直行運(yùn)動(dòng)不明顯;當(dāng)DR較小時(shí)，與直行運(yùn)動(dòng)耦合使粒子表現(xiàn)出一定的活性(同時(shí)具有明顯的直行運(yùn)動(dòng)和轉(zhuǎn)向運(yùn)動(dòng));當(dāng)DR為零時(shí)，活性粒子無轉(zhuǎn)向運(yùn)動(dòng)，只存在直行運(yùn)動(dòng).圖7 給出具有不同轉(zhuǎn)向擴(kuò)散系數(shù)活性粒子組成的活性流體的相對黏度隨剪切速率的變化.活性粒子的體積分?jǐn)?shù)φ=0.5%，活性粒子的直行速度U0=40 μm/s.從圖7 可以看出，轉(zhuǎn)向運(yùn)動(dòng)和直行運(yùn)動(dòng)的耦合作用是導(dǎo)致過渡區(qū)流變曲線非單調(diào)變化的主要原因(結(jié)合圖6 分析)，黏度下降區(qū)活性流體黏度下降的現(xiàn)象是由活性粒子直行運(yùn)動(dòng)所引起的.究其原因是不發(fā)生轉(zhuǎn)向運(yùn)動(dòng)的活性粒子僅具有直行運(yùn)動(dòng)能力，其仍可通過活性粒子間相互碰撞作用形成一定規(guī)模的微集群運(yùn)動(dòng)，造成低剪切速率區(qū)域內(nèi)宏觀黏度的下降;由于形成的微結(jié)構(gòu)穩(wěn)定性較差，剪切速率升高后也容易被破壞，因此在中剪切速率區(qū)域沒有觀察到流變曲線非單調(diào)變化的現(xiàn)象.從圖7 還可以看出，若活性粒子存在轉(zhuǎn)向運(yùn)動(dòng)，低剪切速率區(qū)域活性流體黏度的下降幅度隨DR的增大而減小;在過渡區(qū)，先增稠后稀化的非單調(diào)流變學(xué)現(xiàn)象隨著DR的增大而減弱.這些現(xiàn)象說明若活性粒子的運(yùn)動(dòng)方向改變地過于頻繁，會(huì)對活性流體的非牛頓流體行為有明顯的抑制作用.

圖7 活性粒子轉(zhuǎn)向擴(kuò)散系數(shù)對活性流體相對黏度的影響Fig.7 Effect of rotational diffusion coefficient of active particle on relative viscosity of active fluids

2.2 活性流體特殊流變行為的形成機(jī)理

懸浮液體系的宏觀流變特性是由體系中的介觀結(jié)構(gòu)決定的，通過考察不同剪切速率下活性粒子的聚集行為，可以確定活性流體特殊流變學(xué)行為的形成機(jī)理.活性粒子可自發(fā)形成直行運(yùn)動(dòng)和轉(zhuǎn)向運(yùn)動(dòng)，粒子間相互作用頻繁，容易形成聚集結(jié)構(gòu)，誘導(dǎo)集群運(yùn)動(dòng).圖8 給出活性粒子漲落計(jì)算方法示意圖.如圖8所示，觀察區(qū)被等面積劃分為若干區(qū)域，對每個(gè)區(qū)域內(nèi)粒子數(shù)量計(jì)數(shù)，以每個(gè)區(qū)域內(nèi)粒子漲落(標(biāo)準(zhǔn)差ρ)表征活性粒子的聚集程度.標(biāo)注差ρ的計(jì)算公式如下

圖8 活性粒子漲落計(jì)算方法示意圖Fig.8 Schematic of the calculation method for the fluctuation of active particles

其中，xi為第i個(gè)區(qū)域中的活性粒子數(shù)，為所有區(qū)域中活性粒子數(shù)的平均值，n是劃分區(qū)域的個(gè)數(shù).

圖9 是不同工況下活性流體中粒子漲落隨剪切速率的變化.從圖9 可以看出，活性粒子的漲落現(xiàn)象與活性粒子的體積分?jǐn)?shù)、直行速度和轉(zhuǎn)向擴(kuò)散系數(shù)的大小直接相關(guān);活性粒子體積分?jǐn)?shù)越高、直行速度越快和轉(zhuǎn)向擴(kuò)散系數(shù)越小，活性粒子分布的標(biāo)準(zhǔn)差越大，也就是漲落越顯著.從圖9 還可以看出，隨著外部施加剪切速率的不斷增大，在低剪切速率區(qū)域，粒子分布的標(biāo)準(zhǔn)差緩慢減小;在中剪切速率區(qū)域，粒子分布的標(biāo)準(zhǔn)差急劇的下降;而在高剪切速率區(qū)域，粒子分布的標(biāo)準(zhǔn)差基本不發(fā)生改變.3 個(gè)剪切速率區(qū)域內(nèi)的漲落情況與2.1 節(jié)中活性流體的流變曲線劃分為3 個(gè)階段具有一致性，這說明活性流體特殊的流變行為與活性粒子的漲落密切相關(guān).在低剪切速率區(qū)域，粒子自驅(qū)動(dòng)形成的聚集結(jié)構(gòu)產(chǎn)生的內(nèi)部應(yīng)力可以抵消部分外部剪切應(yīng)力，使得驅(qū)動(dòng)活性流體達(dá)到一定剪切速率的剪切應(yīng)力減小，因此體系的相對黏度ηr＜ 1.隨著外部剪切速率的增大，活性流體中的聚集結(jié)構(gòu)產(chǎn)生的內(nèi)部應(yīng)力相比于外部剪切作用變小，驅(qū)動(dòng)活性流體達(dá)到一定剪切速率所需的剪切力增大，宏觀表現(xiàn)為剪切增稠行為，介觀表現(xiàn)為粒子分布標(biāo)準(zhǔn)差快速減小.在高剪切速率區(qū)域，活性粒子自發(fā)形成的聚集結(jié)構(gòu)會(huì)被逐漸破壞，宏觀表現(xiàn)為剪切稀化行為，介觀表現(xiàn)為粒子分布標(biāo)準(zhǔn)差加速減小.當(dāng)外部施加的剪切作用足夠大時(shí)，活性流體中的聚集結(jié)構(gòu)完全消失，活性粒子的分布趨向于與被動(dòng)粒子類似的隨機(jī)分布，此時(shí)體系可視為一般被動(dòng)流體，具有牛頓流體特性.活性粒子的體積分?jǐn)?shù)、自驅(qū)動(dòng)直行速度和轉(zhuǎn)向擴(kuò)散系數(shù)均通過影響聚集結(jié)構(gòu)的形成進(jìn)而改變體系的宏觀流變特性.

圖9 活性流體中粒子密度標(biāo)準(zhǔn)差與剪切速率的關(guān)系Fig.9 Density standard deviation of particles as function of shear rates in active fluids

3 結(jié)論

本文基于布朗動(dòng)力學(xué)理論，通過擴(kuò)展布朗運(yùn)動(dòng)方程，引入驅(qū)動(dòng)項(xiàng)和隨機(jī)旋轉(zhuǎn)項(xiàng)再現(xiàn)活性粒子的直行和轉(zhuǎn)向運(yùn)動(dòng)，建立可以模擬活性粒子運(yùn)動(dòng)的數(shù)學(xué)模型.采用反向非平衡分子動(dòng)力學(xué)(rNEMD)法計(jì)算活性流體黏度，考察了活性粒子體積分?jǐn)?shù)、直行速度和轉(zhuǎn)向擴(kuò)散系數(shù)對活性流體流變行為的影響規(guī)律，并分析活性粒子的漲落確定活性流體特殊流變行為的形成機(jī)理.結(jié)論如下:

(1)建立的活性流體布朗動(dòng)力學(xué)模型可以很好地描述活性粒子的直行和轉(zhuǎn)向運(yùn)動(dòng)，可以再現(xiàn)活性流體體系中粒子分布不均勻的非平衡特征，rNEMD法可以高效準(zhǔn)確地計(jì)算體系的黏度.

(2)隨著剪切速率的增大，活性流體的流變曲線可劃分為3 個(gè)階段.在低剪切速率區(qū)域?yàn)轲ざ认陆祬^(qū)，活性流體的相對黏度ηr＜ 1;中剪切速率區(qū)域?yàn)檫^渡區(qū)，活性流體的相對黏度呈現(xiàn)先增大后減小的非單調(diào)變化;高剪切速率區(qū)域?yàn)榕ｎD區(qū)，活性流體的黏度與一般被動(dòng)流體類似.活性粒子的體積分?jǐn)?shù)越大、活性粒子的直行速度越快，非單調(diào)流變學(xué)行為越顯著.低剪切速率區(qū)域活性流體黏度下降是因?yàn)榛钚粤Ｗ又毙羞\(yùn)動(dòng)存在引起的，轉(zhuǎn)向運(yùn)動(dòng)的存在是導(dǎo)致中剪切速率區(qū)域流變曲線非單調(diào)變化的主因，活性粒子頻繁改變運(yùn)動(dòng)方向會(huì)抑制活性流體的非牛頓特性.

(3)活性流體的宏觀流變學(xué)特性和粒子的漲落直接相關(guān).活性粒子體積分?jǐn)?shù)越高、直行速度越快和轉(zhuǎn)向擴(kuò)散系數(shù)越小，活性流體中活性粒子越容易產(chǎn)生顯著的漲落;活性粒子的漲落越顯著，低剪切速率區(qū)域的黏度下降現(xiàn)象和中剪切速率區(qū)域的非單調(diào)流變學(xué)行為越明顯.隨著外部剪切流場的增強(qiáng)，活性粒子形成的局部聚集結(jié)構(gòu)會(huì)被外部施加的剪切作用破壞，并在中剪切速率區(qū)域流變曲線出現(xiàn)先增稠后稀化的非單調(diào)行為;當(dāng)剪切速率進(jìn)一步增大時(shí)，活性流體中粒子聚集結(jié)構(gòu)被完全破壞，體系可視為一般被動(dòng)流體.