1、屈服應力流體的粗粒化分子動力學模擬研究摘要：非牛頓屈服應力流體廣泛應用于生活和工業生產中，具有很高的研究價值。利用耗散粒子動力學(DID：Dissipative Particle Dynamics)方法來對屈服應力流體進行模擬研究%首先應用D1D方法構建了一種屈服應力流體模型，并模 擬預測了該流體模型具有Herschel-Bulkley流體的流變特性%通過數值模擬驗證了模型在剪切和伯肅葉流動特性，符合屈 服應力流體的基本特點關鍵詞：非牛頓，屈服應力流體，數值模擬,耗散粒子動力學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)。是一個服從高斯統計分布的隨機波動變量，滿 足=和= (&ik&jl+&il&jk)&(tt)。C，D 和R 分別是這三種力的權重函數，這些權重函數都符合下式的形式：ri在日常生活中，很多物質在通常情況下既不是固體也不是 簡單的牛頓流體，而具有非線性的流變行為，這類物質稱為非牛 頓流體。其中，有一類材料只有在被施加的作用力大于某個閾值 時才會展現出流動的性能，反之則表現出固體的特性，具有這類 臨界流變學特性的材料稱為屈服應力(yield-stress)流體，如: 繪畫顏料、泡沫材料、濕潤的水泥、清潔膏、蛋黃醬和牙膏等)#-2*， 另外在石油工業中，重

5、新啟動凝膠原油管道時能夠準確地估計 最小壓力是十分重要的。由于流動應力閾值的存在，這類材料 具有高度的非線性特點，往往對它流動的動態特性起著極其重 要的作用，具有很高的研究和應用價值。傳統宏觀的數值模擬方法發展起來的理論架構體系對非牛 頓屈服應力問題的預測及分析并不夠準確和細致，為了能從物 質的微觀給出合理的分析驗證，因而需要將微觀方法引入到該 問題中。這些具有獨特的性質，使它們能夠具有聯系宏觀和微觀 的潛力，同時在模擬中既能兼顧計算效率和模型尺度，能夠為微 流體問題的研究工作提供可參考性較高的數值指導。本文采用的是粗?；姆肿觿恿W模擬方法，即：耗散粒子動 力學方法，以實現構建屈服應力流體模

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

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

8、 滿足下式：=) * ，$ =2#kT(5)其中的kpT是系統的Boltzmann溫度，在本文中統一簡稱 為系統的溫度。為了具有高施密特數的流體動力學響應，使得隨 機力與保守力采用相同的權重指數(n=1)0需要指出的是，因為 DPD方法中的隨機力和耗散力滿足式(5)，因此二者構成了 DPD方法內置的恒溫器。本文中模擬用到的核心算法是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）上式中的?。╰

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

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

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

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

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

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