1、屈服應(yīng)力流體的粗粒化分子動力學(xué)模擬研究摘要：非牛頓屈服應(yīng)力流體廣泛應(yīng)用于生活和工業(yè)生產(chǎn)中，具有很高的研究價(jià)值。利用耗散粒子動力學(xué)(DID：Dissipative Particle Dynamics)方法來對屈服應(yīng)力流體進(jìn)行模擬研究%首先應(yīng)用D1D方法構(gòu)建了一種屈服應(yīng)力流體模型，并模 擬預(yù)測了該流體模型具有Herschel-Bulkley流體的流變特性%通過數(shù)值模擬驗(yàn)證了模型在剪切和伯肅葉流動特性，符合屈 服應(yīng)力流體的基本特點(diǎn)關(guān)鍵詞：非牛頓，屈服應(yīng)力流體，數(shù)值模擬,耗散粒子動力學(xué)shear and poiseuille are verified by numerical simulation,a

2、nd itAbstract:Non-Newtonian yield stress fluids are widely used in daily life and industrial production,which shows high research value.This paper uses the Dissipative Particle Dynamics (DPD) method to simulate the yield-stress fluid.A yield-stress fluid model is constructed with DPD method,and the

3、rheological properties of fluid model are consistent with those of Herschel-Bulkley (H-B).The flow characteristics of the model in accorded with the basic characteristics of yield-stress fluid.shear and poiseuille are verified by numerical simulation,and itKeywords:Non-Newtonian,yield-stress fluid,s

4、imulation,DPD(1)(1)(2)。是一個(gè)服從高斯統(tǒng)計(jì)分布的隨機(jī)波動變量，滿 足=和= (&ik&jl+&il&jk)&(tt)。C，D 和R 分別是這三種力的權(quán)重函數(shù)，這些權(quán)重函數(shù)都符合下式的形式：ri在日常生活中，很多物質(zhì)在通常情況下既不是固體也不是 簡單的牛頓流體，而具有非線性的流變行為，這類物質(zhì)稱為非牛 頓流體。其中，有一類材料只有在被施加的作用力大于某個(gè)閾值 時(shí)才會展現(xiàn)出流動的性能，反之則表現(xiàn)出固體的特性，具有這類 臨界流變學(xué)特性的材料稱為屈服應(yīng)力(yield-stress)流體，如: 繪畫顏料、泡沫材料、濕潤的水泥、清潔膏、蛋黃醬和牙膏等)#-2*， 另外在石油工業(yè)中，重

5、新啟動凝膠原油管道時(shí)能夠準(zhǔn)確地估計(jì) 最小壓力是十分重要的。由于流動應(yīng)力閾值的存在，這類材料 具有高度的非線性特點(diǎn)，往往對它流動的動態(tài)特性起著極其重 要的作用，具有很高的研究和應(yīng)用價(jià)值。傳統(tǒng)宏觀的數(shù)值模擬方法發(fā)展起來的理論架構(gòu)體系對非牛 頓屈服應(yīng)力問題的預(yù)測及分析并不夠準(zhǔn)確和細(xì)致，為了能從物 質(zhì)的微觀給出合理的分析驗(yàn)證，因而需要將微觀方法引入到該 問題中。這些具有獨(dú)特的性質(zhì)，使它們能夠具有聯(lián)系宏觀和微觀 的潛力，同時(shí)在模擬中既能兼顧計(jì)算效率和模型尺度，能夠?yàn)槲?流體問題的研究工作提供可參考性較高的數(shù)值指導(dǎo)。本文采用的是粗粒化的分子動力學(xué)模擬方法，即：耗散粒子動 力學(xué)方法，以實(shí)現(xiàn)構(gòu)建屈服應(yīng)力流體模

6、型和完成相關(guān)模擬工作。 1耗散粒子動力學(xué)方法標(biāo)準(zhǔn)的DPD方法首先由Hoogerbrugge和Koelman提 出，是MD自下而上發(fā)展而來的粒子粗粒化的方法之一。對于任 意第i個(gè)DPD粒子符合牛頓第二定律：五=!業(yè)=!dt =!，dt =其中r,、i和！分別為第i個(gè)DPD粒子的位置、速度向量和 受力之和(包括非粒子間的外部作用力)，粒子間作用力方向?yàn)?兩粒子中心連線方向。另外，單個(gè)DPD粒子的質(zhì)量被設(shè)為單位 1。成對的粒子相互間作用力可以分解為三個(gè)在兩粒子中心連線 方向上的分力，即保守力(商)、耗散力(F：)和隨機(jī)力(F：)：=$ (+)以上分力的具體形式如下所示：F= ()?=# ()( !)

7、=$ (WM(3)其中aij,#和$分別是保守力，耗散力和隨機(jī)力的幅值。 /j和Vij分別是第i粒子和第j粒子的相對位移和相對速度。并且 有 / ij = r,3(r“)=從上式(4)中可知三種力都是短程形式的，即在相對距離大 于其對應(yīng)的截?cái)喟霃?/c)時(shí)該部分作用力歸零。式(4)中有兩個(gè) 可選參數(shù)，即截?cái)喟霃?rc)和權(quán)重指數(shù)(n)0 一方面，針對于標(biāo)準(zhǔn) DPD方法，每種力都應(yīng)該給定相應(yīng)的截?cái)嗑嚯x(rc)；另一方面， 根據(jù)Espanol和Warren提出耗散力和隨機(jī)力的權(quán)重函數(shù)及 幅值應(yīng)當(dāng)遵循某些約束以滿足漲落耗散定理。具體來說，如果其 中一個(gè)的系數(shù)和權(quán)重函數(shù)中被任意選擇了，那么另外一個(gè)應(yīng)該

8、 滿足下式：=) * ，$ =2#kT(5)其中的kpT是系統(tǒng)的Boltzmann溫度，在本文中統(tǒng)一簡稱 為系統(tǒng)的溫度。為了具有高施密特?cái)?shù)的流體動力學(xué)響應(yīng)，使得隨 機(jī)力與保守力采用相同的權(quán)重指數(shù)(n=1)0需要指出的是，因?yàn)?DPD方法中的隨機(jī)力和耗散力滿足式(5)，因此二者構(gòu)成了 DPD方法內(nèi)置的恒溫器。本文中模擬用到的核心算法是Velocity-Verlet積分迭代算法：r(t+At) = rO + E (t) + # (At) f(t)!(t+,t)=!(t)+At !(t)!（#+#）=!（&（t+t）,!（t+t）!（t+t）=!（t）+（ #（ !（t）+!（t+t）上式中的！（t

9、）、#（t）和r,（t）分別代表著第i粒子在第t個(gè)時(shí) 刻的加速度矢量，速度矢量和位移矢量;V（t+t）代表第i粒子 在第t+ t時(shí)刻的預(yù)測速度,t取值0.01 &同時(shí),參數(shù)是 DPD模擬方法中引入的經(jīng)驗(yàn)參數(shù)，代表了一些隨機(jī)性質(zhì)的附加 效應(yīng)，在本文中取值0.65。2屈服應(yīng)力流體模型的構(gòu)建與模擬預(yù)測為了構(gòu)建具有屈服應(yīng)力性質(zhì)的非牛頓流體，提出了一種帶 有兩種類型顆粒的DPD流體模型&這里把這兩種顆粒分別簡稱 為Pa和Pb,并且體積分?jǐn)?shù)比例設(shè)置成3:7，如圖1左圖所示&由 于系統(tǒng)初始狀態(tài)的粒子是規(guī)則且分成分布的，所以整體流體系 統(tǒng)并未達(dá)到平衡狀態(tài)，即系統(tǒng)的溫度5bT%1為了能給后期的研 究作好準(zhǔn)備，首

10、先需要將初始狀態(tài)的流體系統(tǒng)模擬到平衡狀態(tài)， 如圖1右圖所示&圖1屈服應(yīng)力DPD流體模型初始構(gòu)建分布和穩(wěn)定狀態(tài)的粒子圖為了預(yù)測流體模型的流變性質(zhì)，需要對此模型施加剪切作 用（取不同的隨機(jī)力參數(shù)#為2.0、3.0、5.0）,使用Irving-Ki rk- wood模型計(jì)算得到應(yīng)力張量，最終可以畫出流體模型的應(yīng)力- 應(yīng)變曲線，如圖2的上圖&圖2 不同#參數(shù)下（M）DPD流體模型的應(yīng)力-應(yīng)變圖從圖2中可以看到該流體存在一個(gè)明顯的屈服應(yīng)力So流變特性是非線性的&我們將此模擬的結(jié)果與Herschel-Bulkley流體9（以下簡稱H-B流體）理論模型進(jìn)行對比驗(yàn)證， 得到了很好的吻合度，如圖2的下圖所示&另

11、外,H-B流體的流變特性可以用下式表達(dá)：T = T0+k%（7）式中的指數(shù)n的取值范圍為0.20.8,具有剪切變稀的特 點(diǎn)&從圖2中我們可以看到在不同隨機(jī)力參數(shù)#情況下，擬合曲 線的n的取值都在以上的范圍內(nèi)，分別為0.31 $0.24和0.45。 3數(shù)值模擬結(jié)果 3.1剪切流驗(yàn)證本節(jié)中我們將對模型進(jìn)行剪切作用流動的模擬驗(yàn)證&我們 將圖3中的圖像根據(jù)模擬設(shè)置分成靜置形成、施加剪切破壞和 靜置恢復(fù)結(jié)構(gòu)等三個(gè)階段，以便于我們做進(jìn)一步的闡述分析&在 靜置形成結(jié)構(gòu)階段，流體粒子經(jīng)過充分地混合后達(dá)到系統(tǒng)平衡 狀態(tài),形成取向比較一致的拓?fù)浣Y(jié)構(gòu)，如圖3a,在粒子Pb的體 積分?jǐn)?shù)分布云圖3A中也可以明顯地看出，

12、各區(qū)域體積分?jǐn)?shù)的 差異較大，且呈現(xiàn)一定的規(guī)律性&在施加剪切階段，需要對施加 大小不同的剪切速率下的工況分別進(jìn)行分析，例如圖3B和b 是在施加剪切速率=0.5情況下得到的結(jié)果，可以看到各粒子 的分布結(jié)構(gòu)在靠近壁面處先被剪切流動所破壞；當(dāng)施加的剪切 速率增大到1.0的時(shí)候，各粒子的分布結(jié)構(gòu)則被完全遭到破壞, 粒子和體積分?jǐn)?shù)趨向于均勻分布，如圖3C和c所示&在靜置恢 復(fù)結(jié)構(gòu)階段，流體系統(tǒng)未受到任何的剪切作用，粒子又重新成恢 復(fù)原有的拓?fù)浣Y(jié)構(gòu)狀態(tài)，體積分?jǐn)?shù)分布亦然，這在圖3D和d中 得到了充分地展示&且圖4剪切作用力流場的速度剖面圖：（左）屈服應(yīng)力流體（右）牛頓流體圖3 （M）DPD模擬中流場的體積分

13、數(shù)云圖和粒子結(jié)構(gòu)圖的變化情況 根據(jù)剪切流場的速度剖面圖，非牛頓屈服應(yīng)力流體隨著剪 切速率增大的同時(shí)速度剖面形狀會出現(xiàn)一個(gè)明顯的變化，具體 來說，在小于轉(zhuǎn)折剪切速率的情況下，速度剖面近似地呈現(xiàn)一種 波浪號的特點(diǎn)，即在靠近中心線附近的粒子層之間沒有出現(xiàn)剪 切流動特點(diǎn)，出現(xiàn)“塊狀”的流動；而在大于轉(zhuǎn)折剪切速率的條件 下,其速度剖面會出現(xiàn)類似于圖4（右）中的牛頓流體模型那樣 的線性剪切流動現(xiàn)象&這種特殊情況的出現(xiàn)也可以利用屈服應(yīng) 力的存在來解釋，當(dāng)剪切速率大于某一特定值時(shí)，流體各部分的 應(yīng)力值大于屈服應(yīng)力且圖4剪切作用力流場的速度剖面圖：（左）屈服應(yīng)力流體（右）牛頓流體3.2泊肅葉流驗(yàn)證對屈服應(yīng)力流體

14、施加多種大小不同的體積力，用以研究流 體的泊肅葉流動情況$待流體模型到達(dá)穩(wěn)定狀態(tài)以后，對兩種類 型的流體顆粒在x軸方向上施加一個(gè)額外恒定的體積力，即 g=常數(shù)$另外，將x軸和y軸方向的邊界條件設(shè)置成周期性邊 界條件。經(jīng)過一段時(shí)間的模擬后，可以繪制流體在x軸方向上的 速度分量隨著z軸坐標(biāo)的變化曲線，得到圖5（上）$通過分析泊肅葉流動的速度剖面，我們可以很容易地得到速 度剖面可以分成兩個(gè)部分，即拋物線形式流動部分和平臺形式流 動部分$為了更加直觀地展示這兩種流動區(qū)域，我們繪制得到施 加的體積力為1.00（即gx=1.00）的流動模型的速度云圖，如圖5 （下）所示。上述的平臺流動部分位于整個(gè)流域的中心線附近，且 對稱分布，而拋物線形式流動部分則位于靠近壁面的區(qū)域$ 的區(qū)域大小比例也隨著體積力的變化而變化。具體來說，施加的 體積力越大，