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

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

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

過(guò)往的模擬研究為了簡(jiǎn)化計(jì)算，通常會(huì)忽略一些影響因素，例如粒子間作用力、活性粒子的轉(zhuǎn)向運(yùn)動(dòng)和直行速度的變化等，因此會(huì)出現(xiàn)和實(shí)驗(yàn)結(jié)果不一致的情況.此外，還存在不能發(fā)揮模擬計(jì)算單因素分析優(yōu)勢(shì)的問(wèn)題，僅考慮活性粒子體積分?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ù)對(duì)活性流體流變行為的影響規(guī)律，并且分析活性粒子的漲落規(guī)律，探究活性流體特殊流變行為的形成機(jī)理.

1 模型的建立和求解

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

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

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

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

活性粒子為亞微米或微米尺寸，其運(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)能.因此，除了流體分子對(duì)粒子的碰撞和摩擦作用以外，活性粒子還在自驅(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來(lái)描述.

圖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)力張量的分量.然而，對(duì)于活性流體這樣的非平衡體系，使用NEMD 法計(jì)算會(huì)得到更加準(zhǔn)確的結(jié)果.NEMD 法的思路是對(duì)流體施加一個(gè)剪切作用，計(jì)算流體在剪切作用下的響應(yīng)，借助計(jì)算動(dòng)量通量和速度梯度的比值確定流體的黏度[29-30].Müller-Plathe[31]對(duì)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ā)生了交換.隨著模擬過(guò)程的進(jìn)行，可以在體系中產(chǎn)生一個(gè)剪切速率曲線來(lái)計(jì)算流體的黏度.該方法也被稱為反向非平衡分子動(dòng)力學(xué)(reverse non-equilibrium molecular dynamics，rNEMD)法，這是因?yàn)镹EMD 法是通過(guò)對(duì)體系施加剪切作用，然后對(duì)應(yīng)力張量的非對(duì)角元(正比于動(dòng)量通量)進(jìn)行統(tǒng)計(jì)計(jì)算，而rNEMD 法是先施加動(dòng)量通量，剪切速率曲線是流體響應(yīng)的結(jié)果.

在布朗動(dòng)力學(xué)的模擬過(guò)程中，通過(guò)對(duì)流體速度場(chǎng)、特征剪切應(yīng)力和特征剪切速率等物理量的統(tǒng)計(jì)計(jì)算，就可以獲得流體的黏度值.為驗(yàn)證黏度計(jì)算方法的可靠性，采用簡(jiǎn)單LJ 流體進(jìn)行測(cè)試，其理論黏度η*=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，更接近該流體的理論黏度值，且相對(duì)于其他兩種方法計(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 算法，粒子間作用勢(shì)能選擇LJ 勢(shì)函數(shù)，運(yùn)算時(shí)的物理量使用LJ 單位，邊界處采用周期性邊界條件.本文為了考察活性粒子體積分?jǐn)?shù)、直行速度和轉(zhuǎn)向擴(kuò)散系數(shù)對(duì)活性流體流變行為的影響規(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 (無(wú)直行運(yùn)動(dòng))、20 μm/s，30 μm/s 和40 μm/s，活性粒子的轉(zhuǎn)向行為由轉(zhuǎn)向擴(kuò)散系數(shù)DR反映，大小分別為0 s-1(無(wú)轉(zhuǎn)向運(yùn)動(dòng))、1.0 × 10-4s-1、5.0 × 10-4s-1和1.0 × 10-3s-1.

通過(guò)求解描述活性粒子運(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í)間步長(zhǎng)Δt=0.001τ.使用大量獨(dú)立軌跡的時(shí)間平均來(lái)獲得系統(tǒng)平均，這對(duì)于體系宏觀物理量的計(jì)算是必要的.

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

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

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

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

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

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

圖5 給出不同體積分?jǐn)?shù)活性流體的相對(duì)黏度ηr隨剪切速率的變化.活性粒子的自驅(qū)動(dòng)速度U0=40 μm/s，轉(zhuǎn)向擴(kuò)散系數(shù)DR=0.000 1 s-1.從圖5 可以看出，對(duì)于體積分?jǐn)?shù)為0.05%的活性流體，相對(duì)黏度隨剪切速率增大沒(méi)有明顯改變;而體積分?jǐn)?shù)φ≥0.1%的活性流體呈現(xiàn)非牛頓流體特性，流變曲線具有非單調(diào)性，可以劃分為3 個(gè)區(qū)域，具體包括:(1) 黏度下降區(qū)(＜～ 10 s-1)，低剪切速率區(qū)域內(nèi)活性流體的相對(duì)黏度ηr＜1，即活性流體的黏度低于未加入活性粒子的流體黏度，并且活性粒子的體積分?jǐn)?shù)越高，黏度下降越明顯;(2) 過(guò)渡區(qū)(～ 10 s-1＜＜～ 50 s-1)，中剪切速率區(qū)域內(nèi)活性流體的黏度出現(xiàn)先增大后減小的變化趨勢(shì)，并且存在一個(gè)臨界剪切速率c，此時(shí)活性流體的黏度最大，臨界剪切速率c出現(xiàn)在剪切速率=20 s-1附近;(3) 牛頓區(qū)(＞～ 50 s-1)，高剪切速率區(qū)域內(nèi)活性流體的相對(duì)黏度η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)黏度變化的過(guò)渡區(qū)域，黏度值小幅度增大.

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

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

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

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

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

懸浮液體系的宏觀流變特性是由體系中的介觀結(jié)構(gòu)決定的，通過(guò)考察不同剪切速率下活性粒子的聚集行為，可以確定活性流體特殊流變學(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ū)域，對(duì)每個(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è)階段具有一致性，這說(shuō)明活性流體特殊的流變行為與活性粒子的漲落密切相關(guān).在低剪切速率區(qū)域，粒子自驅(qū)動(dòng)形成的聚集結(jié)構(gòu)產(chǎn)生的內(nèi)部應(yīng)力可以抵消部分外部剪切應(yīng)力，使得驅(qū)動(dòng)活性流體達(dá)到一定剪切速率的剪切應(yīng)力減小，因此體系的相對(duì)黏度η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ù)均通過(guò)影響聚集結(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é)理論，通過(guò)擴(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ù)對(duì)活性流體流變行為的影響規(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ū)，活性流體的相對(duì)黏度ηr＜ 1;中剪切速率區(qū)域?yàn)檫^(guò)渡區(qū)，活性流體的相對(duì)黏度呈現(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é)行為越明顯.隨著外部剪切流場(chǎng)的增強(qiáng)，活性粒子形成的局部聚集結(jié)構(gòu)會(huì)被外部施加的剪切作用破壞，并在中剪切速率區(qū)域流變曲線出現(xiàn)先增稠后稀化的非單調(diào)行為;當(dāng)剪切速率進(jìn)一步增大時(shí)，活性流體中粒子聚集結(jié)構(gòu)被完全破壞，體系可視為一般被動(dòng)流體.