耗散粒子動(dòng)力學(xué)中Lees－Edwards邊界條件的實(shí)施

陳 碩，金亞斌，張明焜，尚 智

（1.同濟(jì)大學(xué) 航空航天與力學(xué)學(xué)院，上海200092；2.高性能計(jì)算研究所，A＊STAR，新加坡138632）

Lees－Edwards邊界條件（簡(jiǎn)稱(chēng)LE）最初在分子動(dòng)力學(xué)中模擬簡(jiǎn)單剪切流動(dòng)［1］，該邊界條件的特點(diǎn)是在不需要固體壁面拖動(dòng)的情況下在準(zhǔn)無(wú)限大流體系統(tǒng)中獲得定常的剪切率.

基于粒子的數(shù)值模擬方法中，密集排列的粒子一般會(huì)對(duì)近壁面區(qū)域內(nèi)流體密度、溫度等參量有一定的影響，進(jìn)而影響到流動(dòng)的主體.LE 方法在邊界附近的區(qū)域內(nèi)則不會(huì)出現(xiàn)物理參數(shù)的數(shù)值反常波動(dòng).LE 邊界條件其原理也適合于基于粒子的其他模擬方法，最近Lorenz等［2］，Wagner等［3］研究了LE邊界條件在格子－波爾茲曼方法（簡(jiǎn)稱(chēng)LB）中的實(shí)施，并模擬了懸浮液的剪切流動(dòng)等.

將從微觀(guān)、介觀(guān)到宏觀(guān)不同層次及尺度的模擬聯(lián)系起來(lái)正成為計(jì)算機(jī)模擬領(lǐng)域最富挑戰(zhàn)性的重要課題.介觀(guān)模擬更是一項(xiàng)為實(shí)現(xiàn)“第一原理預(yù)計(jì)物質(zhì)宏觀(guān)性質(zhì)”奠定物理基礎(chǔ)的關(guān)鍵技術(shù)［4］，而發(fā)展適宜的介觀(guān)尺度的流體模擬技術(shù)是近年來(lái)國(guó)際計(jì)算機(jī)模擬領(lǐng)域的熱點(diǎn)［5］，其特征為在分子尺度上進(jìn)行粗?；幚恚春雎宰罴?xì)節(jié)的分子結(jié)構(gòu)），而在介觀(guān)尺度上捕捉流體的物理特性.

耗散粒子動(dòng)力學(xué)（簡(jiǎn)稱(chēng)DPD）是一門(mén)新興的基于粒子的介觀(guān)尺度數(shù)值模擬技術(shù)，于1992 年首次由Hoogerbrugge和Koelman提出［6］.盡管DPD 不如LB計(jì)算省時(shí)，但是DPD 比LB 靈活得多，尤其不會(huì)象LB那樣在很多情況下出現(xiàn)數(shù)值不穩(wěn)定性［7］.另一方面，DPD 又可以理解為宏觀(guān)的微分流動(dòng)控制方程在小尺度上的隨機(jī)描述，從這個(gè)意義上說(shuō)，DPD 方法是連接微觀(guān)分子動(dòng)力學(xué)方法和宏觀(guān)流體力學(xué)方法的一座橋梁.

DPD 已經(jīng)成功地應(yīng)用于對(duì)復(fù)雜流體和復(fù)雜流動(dòng)過(guò)程的模擬［8－13］，其中包括早期的DPD算法結(jié) 合LE邊界條件獲得簡(jiǎn)單剪切流動(dòng)［6，14］.隨著DPD 中計(jì)算方法的改進(jìn)，在DPD 中實(shí)施LE邊界條件，需要將LE條件與DPD 中所采用的velocity－Verlet算法以及元胞分割方法等結(jié)合起來(lái)考慮.本文結(jié)合上述算法詳細(xì)討論了LE條件在DPD 方法中的實(shí)施過(guò)程，并對(duì)在LE 邊界條件下所獲得的速度、密度、壓強(qiáng)、剪切應(yīng)力等參數(shù)的分布進(jìn)行了比較分析與討論.通過(guò)提高耗散力系數(shù)來(lái)模擬較高的Schmidt數(shù)流體，結(jié)果表明在高耗散率條件下，在耗散粒子動(dòng)力學(xué)中應(yīng)用Lees－Edwards邊界條件仍然有效.

1 耗散粒子動(dòng)力學(xué)

1.1 耗散粒子動(dòng)力學(xué)理論模型

基于牛頓運(yùn)動(dòng)方程，DPD 系統(tǒng)中粒子的位置和速度隨時(shí)間的變化規(guī)律為

式中：ri和vi為第i個(gè)粒子的位置和速度矢量；t為時(shí)間；fi為所有其他粒子（不包括粒子i本身）作用在粒子i上的力；Fe為粒子所受到的其他外力如高分子鏈的彈性力等.DPD 系統(tǒng)中單位質(zhì)量為單個(gè)粒子的質(zhì)量，因此作用在每個(gè)粒子上的力的大小等于粒子的加速度值.粒子之間的動(dòng)力相互作用由耗散力和隨機(jī)力組成.

fi包括保守力FC，ij、耗散力FD，ij和隨機(jī)力FR，ij三部分，其中每一部分都表示粒子之間的兩兩相互作用.

式中的求和符號(hào)對(duì)一定的截?cái)喟霃絩c內(nèi)所有的不包括粒子i的其他粒子進(jìn)行求和.如果粒子之間的距離rij＞rc，則粒子之間的相互作用為零.

保守力FC，ij為作用在2個(gè)粒子質(zhì)心連線(xiàn)方向上的粒子之間的排斥力，可表示為

式中：aij為保守力系數(shù)，表示粒子i與粒子j間最大的排斥力值；rij＝ri－rj，rij＝｜rij｜，r＾ij＝rij／｜rij｜.粒子j作用在粒子i上的耗散力（或者阻力）FD，ij可以表達(dá)為

式中：γ為耗散系數(shù)；wD為與粒子間距離rij有關(guān)的權(quán)函數(shù)，當(dāng)rij≥rc時(shí)，wD＝0；vij為粒子之間相對(duì)速度，vij＝vi－vj.隨機(jī)力FR，ij可表達(dá)為

式中：σ為隨機(jī)力系數(shù)；wR為與粒子之間距離rij有關(guān)的權(quán)函數(shù)，稱(chēng)為隨機(jī)力權(quán)函數(shù)，當(dāng)rij≥rc時(shí)，wR＝0；ξij為 隨 機(jī) 變 量，其 平 均 值 為 零 且 方 差 為Δt－1，Δt為時(shí)間步長(zhǎng).與保守力一樣，隨機(jī)力和耗散力也作用在粒子質(zhì)心的連線(xiàn)方向.

Espa?ol等［15］研究表明，式（4）和（5）中的2 個(gè)權(quán)函數(shù)相互關(guān)聯(lián)，兩者只能任取一個(gè)，其相互關(guān)系為

γ和σ也相互關(guān)聯(lián)，兩者也只能任取一個(gè)

式中：kB為Boltzmann常 數(shù)，T為 溫 度，kBT為 系 統(tǒng)的Boltzmann溫度，與系統(tǒng)的漲落－耗散理論相類(lèi)似［16］.將kBT作為能量的單位，則有

Fan等［17］的研究表明，采用平方根的權(quán)函數(shù)有利于提高DPD 系統(tǒng)的Schmidt數(shù)，因此本文中采用改進(jìn)的權(quán)函數(shù)表達(dá)式，即rc取單位長(zhǎng)度，即rc＝1.

1.2 粒子間作用力的計(jì)算

在DPD 中如何恰當(dāng)選取方法計(jì)算粒子間相互作用力將在很大程度上影響整個(gè)模擬所花費(fèi)的時(shí)間，元胞分割法（cell subdivision）［18］不 失為一種高效的方法.如圖1所示，計(jì)算區(qū)域被劃分為一系列的元胞區(qū)域，每個(gè)元胞的邊長(zhǎng)均大于rc.計(jì)算區(qū)域內(nèi)的各DPD 粒子處于不同的元胞區(qū)域內(nèi).顯然，粒子之間的相互作用僅僅發(fā)生在處于同一個(gè)元胞內(nèi)的粒子之間以及分別處于2個(gè)緊鄰的元胞內(nèi)的粒子之間；否則，粒子之間的相互距離必定大于rc，從而可以不考慮其相互作用.因此，在計(jì)算某個(gè)粒子與周?chē)Ｗ拥南嗷プ饔脮r(shí)，可以只計(jì)算該粒子與其屬于同一個(gè)元胞內(nèi)的其他粒子之間的相互作用以及該粒子與緊鄰的周?chē)麅?nèi)的粒子的相互作用.這樣的方法，避免了所求解的粒子與計(jì)算區(qū)域內(nèi)的粒子之間全部求解相互作用力，可以大幅節(jié)約計(jì)算工作量.

由于粒子間相互作用力的對(duì)稱(chēng)性，因此可以?xún)H考慮相鄰的元胞中的一半數(shù)目，對(duì)于二維問(wèn)題來(lái)說(shuō)，是5個(gè)元胞（含中心元胞本身）；對(duì)于三維問(wèn)題來(lái)說(shuō)，是14個(gè)元胞（同樣也含中心元胞本身）.對(duì)于周期性邊界條件，相鄰的鏡像元胞也必須計(jì)算在內(nèi).這些細(xì)節(jié)，在與LE邊界條件的結(jié)合中必須仔細(xì)考慮.

1.3 數(shù)值積分算法

圖1 元胞分割法示意Fig.1 Schematic diagram of cell subdivision method

目前DPD 中已采用改進(jìn)的velocity－Verlet算法［19］來(lái)更新每個(gè)粒子在新時(shí)刻（即t＋Δt時(shí)刻）的位置、速度以及受力情況.如下所示：ri（t＋Δt）＝式中：為粒子i在t＋Δt時(shí)刻的預(yù)估速度；λc為經(jīng)驗(yàn)系數(shù).如果粒子的受力情況與粒子間的相對(duì)速度無(wú)關(guān)，則λc＝0.5.Groot和Warren［19］發(fā)現(xiàn)，DPD方法中λc的最優(yōu)值為0.65.當(dāng)采用該最優(yōu)值同時(shí)取ρ＝3，σ＝3時(shí)，對(duì)平衡系統(tǒng)的模擬中計(jì)算的時(shí)間步長(zhǎng)可以增加到Δt＝0.06而不失去對(duì)溫度的控制.

2 LE邊界條件與DPD中算法的結(jié)合

2.1 LE邊界條件

LE 邊界條件可模擬周期性的簡(jiǎn)單剪切流動(dòng)［20］.圖2為周期性系統(tǒng)中簡(jiǎn)單剪切流動(dòng)示意圖，圖中L為周期性區(qū)域在z方向的長(zhǎng)度，ux為x方向的流體剪切速度.假定所研究的區(qū)域內(nèi)僅有2個(gè)粒子.在粒子運(yùn)動(dòng)的過(guò)程中，受到處于本區(qū)域內(nèi)其他粒子及周期性鏡像區(qū)域內(nèi)粒子的作用.由于剪切率Dx＝?ux／?z，鏡像區(qū)域各自角點(diǎn)的x方向速度u（r，t）與z成正比，即u（r，t）＝iDxz，其中i為x方向單位矢量.

圖2 LE 邊界條件示意Fig.2 Schematic diagram of LE boundary condition

在粒子運(yùn)動(dòng)過(guò)程中，如果粒子跑出了計(jì)算區(qū)域，該粒子將由其周期性鏡像粒子代替而重新進(jìn)入該計(jì)算區(qū)域.然而在剪切流動(dòng)條件下，粒子穿過(guò)z＝0面或者z＝L面時(shí)，其對(duì)應(yīng)的周期性鏡像粒子并不會(huì)有同樣的流體速度以及同樣的x坐標(biāo).跑出計(jì)算區(qū)域的粒子以周期性鏡像的身份再進(jìn)入該區(qū)域時(shí)如果賦予特定的速度和x坐標(biāo)可產(chǎn)生剪切流動(dòng).整個(gè)流動(dòng)過(guò)程中其空間上是均勻的，不會(huì)感受到邊界的存在，這是LE邊界條件的優(yōu)點(diǎn).

具體來(lái)說(shuō)，在模擬區(qū)域中，令其角點(diǎn)O的流體速度為零.所模擬區(qū)域中DPD 粒子的速度由兩部分構(gòu)成，即DPD 粒子的熱速度ci和流體速度u（ri），即

簡(jiǎn)單剪切中由于流體速度是只以z為自變量的函數(shù)，這里只需要考慮在z方向上粒子穿越邊界的情況.零時(shí)刻ri的鏡像矢量r′i在ri＋kL處，ri的鏡像矢量r″i在ri－kL處.t時(shí)間后，該粒子及鏡像的位置分別為其中，ci和zi（包括其鏡像）都是時(shí)間的函數(shù).由于粒子的熱速度及其所有鏡像粒子的熱速度都是相同的，因此，那么同理kL－iDxLt.如果ri（t）跑出了下邊界，將被其鏡像代替

如果ri（t）跑出了上邊界，將被其鏡像r″i（t）代替

2.2 DPD中的算法與LE邊界條件的結(jié)合

如圖3所示，模擬區(qū)域顯示x，z方向，y方向垂直紙面向里，圖中的方格代表元胞.由于流動(dòng)存在剪切率Dx，所以粒子離開(kāi)下邊界而以鏡像粒子的身份從上邊界被引入時(shí)x方向的速度要加上DxLz，即速率變?yōu)閡x＋DxLz.由于x方向速度的差異，一定時(shí)間后在上邊界被引入的x方向的位置也會(huì)相應(yīng)改變?chǔ)，也就是圖形上部實(shí)線(xiàn)粒子與虛線(xiàn)粒子x方向距離差Δx＝ntVxΔt－［ntVxΔt／Lx］Lx，其中nt為時(shí)間步數(shù)，Δt為時(shí)間，［］為取整符號(hào)，Lx為模擬區(qū)域x方向的總長(zhǎng)度.從計(jì)算區(qū)域下部跑出去的粒子將從上邊界處一個(gè)新的位置進(jìn)入計(jì)算區(qū)域，從而能夠保證能量守恒［1］.相應(yīng)地在計(jì)算粒子間相互作用時(shí)其相鄰的粒子也需要更新，因此DPD 所采用的元胞分割法中的所有細(xì)節(jié)也必須與此相對(duì)應(yīng).

圖3 DPD 方法中LE 邊界條件的實(shí)施Fig.3 Implementation of LE boundary condition in DPD

粒子從計(jì)算區(qū)域底部離開(kāi)的算法如下：r′x＝必須強(qiáng)調(diào)的是改進(jìn)的velocity－Verlet算法引入一個(gè)中間速度，對(duì)于粒子從計(jì)算

3 結(jié)果及討論

3.1 LE邊界條件下DPD簡(jiǎn)單剪切流動(dòng)模擬

首先，利用上述DPD 與LE相結(jié)合的方法，模擬了大小為－15.000≤x＜15.000，－5.000≤y＜5.000，－11.125≤z＜11.125區(qū)域內(nèi)的簡(jiǎn)單剪切流動(dòng).總共采用了27 600個(gè)DPD 粒子，由此對(duì)應(yīng)的流體密度為4.13.在開(kāi)始的2 000時(shí)間步長(zhǎng)內(nèi)，令流體粒子自由運(yùn)動(dòng)，不施加LE 邊界條件，直到達(dá)到平衡態(tài).從2 000時(shí)間步開(kāi)始，對(duì)系統(tǒng)施加LE 邊界條件.時(shí)間步長(zhǎng)設(shè)置為0.02.在z方向設(shè)置40 個(gè)統(tǒng)計(jì)格子，所有的局部參數(shù)將由對(duì)各個(gè)格子內(nèi)的數(shù)據(jù)樣本進(jìn)行10 000時(shí)間步長(zhǎng)的平均而獲得.經(jīng)過(guò)45 050步的計(jì)算模擬結(jié)果如圖4所示.

圖4 所示為2 種剪切率條件下的簡(jiǎn)單剪切流動(dòng)，剪切率分別為0.2和0.4.所得簡(jiǎn)單剪切流動(dòng)速度的線(xiàn)性化程度很好，更為重要的是，在邊界區(qū)域，也即在z＝11.125和z＝－11.125附近區(qū)域，模擬的速度沒(méi)有反常波動(dòng)，即感受不到邊界的存在，說(shuō)明模擬結(jié)果與預(yù)想的一致.而其各自的斜率剛好等于剪切率，從而驗(yàn)證了在DPD 中實(shí)施LE 邊界條件的正確性.相應(yīng)地，還統(tǒng)計(jì)了系統(tǒng)內(nèi)密度和溫度分布.圖5所示在剪切率為0.2時(shí)的密度和溫度分布，顯示在邊界附近兩者均無(wú)反常，整個(gè)系統(tǒng)的密度保持在4.1附近，而溫度也控制在1.0附近.

在粒子方法中，如果采用固體壁面模型，在實(shí)施無(wú)滑移和不可穿透固體壁面條件的同時(shí)，很可能會(huì)引入壁面效應(yīng)，對(duì)近壁面區(qū)域流體參量產(chǎn)生一定的影響.Fan等［17］提出的壁面模型即采用了“凍住”的粒子代表壁面，在壁面附近的流體區(qū)域設(shè)置一薄層，在這一薄層內(nèi)流體流動(dòng)保持無(wú)滑移條件，并令薄層內(nèi)的粒子速度方向隨機(jī)分布，速度的大小與給定的壁面溫度相對(duì)應(yīng)，但各個(gè)粒子速度的矢量平均值為零（與靜止的壁面相對(duì)應(yīng)）.進(jìn)一步，令在該薄層內(nèi)的粒子以一定速度離開(kāi)壁面，從而保證粒子不穿透壁面.其離開(kāi)速度式中：vR為隨機(jī)速度矢量；n為壁面單位法向矢量（指向流體），詳見(jiàn)文獻(xiàn)［17］.

圖6為DPD 方法中采用Fan等［17］提出的壁面模型模擬剪切率為0.2時(shí)系統(tǒng)的有關(guān)物理量分布與LE邊界條件所得結(jié)果的對(duì)比，圖6a顯示在壁面附近固體壁面模型產(chǎn)生了一定的密度波動(dòng).DPD 中壓力和剪切應(yīng)力的計(jì)算可參閱文獻(xiàn)［10］.由圖6b～6d可見(jiàn)，對(duì)LE邊界條件來(lái)說(shuō)，本文模擬得到的DPD 系統(tǒng)的溫度、壓力和剪切應(yīng)力在邊界附近過(guò)渡平穩(wěn)，而采用固體壁面模型所得結(jié)果則有較大的波動(dòng).

3.2 應(yīng)用LE邊界條件模擬較高施密特?cái)?shù)流體的剪切流動(dòng)

Schmidt數(shù)Sc是動(dòng)量擴(kuò)散與流體自擴(kuò)散的比值，本文中采用改進(jìn)的權(quán)函數(shù)表達(dá)式（式（9）），Sc數(shù)可表示為［12］：Sc＝0.5＋（2πγρr4c）2（1999kBT），式中：γ為耗散系數(shù)，ρ為密度.傳統(tǒng)的DPD 系統(tǒng)所描述的流體Sc數(shù)偏低，提高耗散系數(shù)是提高流體的Schmidt數(shù)的一種有效途徑，從而使所模擬的流體更接近實(shí)際流體［12］.本文中模擬了σ＝14.14，即γ＝100情況下，時(shí)間步長(zhǎng)為0.02、計(jì)算區(qū)域?yàn)?0.00×10.00×22.25時(shí)系統(tǒng)的速度分布.如圖7所示，在較高施密特?cái)?shù)條件下，剪切率為0.2時(shí)，系統(tǒng)內(nèi)的速度分布仍然是線(xiàn)性的，流動(dòng)沒(méi)有邊界效應(yīng)，說(shuō)明速度梯度在整個(gè)空間區(qū)域中是均勻的.

圖7 γ＝100時(shí)溫度及速度分布Fig.7 Temperature and velocity profile whenγ ＝100

文獻(xiàn)［21］中顯示，當(dāng)γ＝100時(shí)，所模擬的剪切流動(dòng)在整個(gè)系統(tǒng)中不是均勻分布，如圖8所示，其線(xiàn)性速度分布出現(xiàn)了轉(zhuǎn)折.文獻(xiàn)［21］認(rèn)為在高的耗散率條件下，LE 邊界條件會(huì)帶來(lái)嚴(yán)重的誤差，文獻(xiàn)［21］將此誤差歸結(jié)為DPD 中耗散力的作用，并需要進(jìn)行修正才能獲得線(xiàn)性的速度分布.不過(guò)筆者認(rèn)為，LE邊界條件模擬了一個(gè)準(zhǔn)無(wú)限大區(qū)域內(nèi)的剪切流動(dòng)，其算法固有的特點(diǎn)是感受不到邊界的存在.因此，DPD 中耗散力的影響在整個(gè)區(qū)域內(nèi)應(yīng)當(dāng)均勻的，而速度分布也應(yīng)一直保持線(xiàn)性.圖7 的結(jié)果也支持這個(gè)觀(guān)點(diǎn).因此本文作者認(rèn)為，高耗散率條件下，在耗散粒子動(dòng)力學(xué)中應(yīng)用Lees－Edwards邊界條件仍然有效.

圖8 文獻(xiàn)［21］中γ＝100時(shí)速度分布Fig.8 Velocity profile in Ref.［21］whenγ＝100

Groot和Warren［19］指出隨著σ的增加，可能會(huì)引起系統(tǒng)溫度的不穩(wěn)定，因此有必要檢驗(yàn)σ＝14.14也即γ＝100情況下系統(tǒng)的溫度，模擬結(jié)果如圖7，此時(shí)系統(tǒng)溫度依然能夠很好地控制在0.9～1.0之間.

4 結(jié)論

詳細(xì)討論了DPD 方法中實(shí)施LE邊界條件的過(guò)程，涉及DPD 方法中改進(jìn)的velocity－Verlet算法、元胞分割法等與LE 邊界條件的結(jié)合，并將DPD 方法中現(xiàn)有的一種固體壁面模型與LE 條件所模擬的結(jié)果進(jìn)行了比較.LE條件模擬所得到的速度、密度、溫度、壓強(qiáng)以及應(yīng)力分布等均為線(xiàn)性分布，說(shuō)明流體系統(tǒng)中各參量是均勻的，感受不到計(jì)算邊界的存在，符合預(yù)期.在模擬高耗散率流體的剪切流動(dòng)時(shí)，所得到的速度分布仍然為線(xiàn)性分布，表明LE 邊界條件在高耗散率流體的模擬中仍然適用，得出了與文獻(xiàn)［21］不同的觀(guān)點(diǎn).LE 條件在DPD 方法中的實(shí)施可為復(fù)雜流體的剪切流動(dòng)研究提供便利.

［1］ Lees A W，Edwards S F.The computer study of transport processes under extreme conditions［J］.J Phys C：Solid State Phys，1972，5：1921.

［2］ Lorenz E，Hoekstra A G.Lees－Edwards boundary conditions for lattice Boltzmann suspension simulations ［J］.Physical Review E，2009，79：036706.

［3］ Wagner A J，Pagonabarraga I.Lees－ Edwards boundary conditions for Lattice Boltzmann ［J］.Journal of Statistical Physics，2002，107：521.

［4］ 李紅霞，強(qiáng)洪夫.耗散粒子動(dòng)力學(xué)模擬方法的發(fā)展和應(yīng)用［J］.力學(xué)進(jìn)展，2009，39（2）：165.LI Hongxia，QIANG Hongfu.Development and application of the simulation method of dissipative particle dynamics ［J］.Advances in Mechanics，2009，39（2）：165.

［5］ 梁文平，楊俊林，陳擁軍，等.新世紀(jì)的物理化學(xué)：學(xué)科前沿與展望［M］.北京：科學(xué)出版社，2004.LIANG Wenping，YANG Junlin，CHEN Yongjun，et al.The new century physical chemistry：subject front and prospects［M］.Beijing：Science Press，2004.

［6］ Hoogerbrugge P J，Koelman J M V A.Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics［J］.Europhysics Letters，1992，19（3）：155.

［7］ Liu M B，Meakin P，Huang H.Dissipative particle dynamics with attractive and repulsive particle－particle interactions［J］.Physics of Fluids，2006，18：017101.

［8］ Yuan S L，Cai Z T，Xu G Y.Mesoscopic simulation of aggregates in surfactant／oil／water systems ［J］.Chinese Journal of Chemistry，2003，21（2）：112.

［9］ Kong Y，Manke C W，Madden W G，et al.Effect of solvent quality on the conformation and relaxation of polymers via dissipative particle dynamics［J］.Journal of Chemical Physics，1997，107：592.

［10］ Chen S，Phan－Thien N，F(xiàn)an X J，et al.Dissipative particle dynamics simulation of polymer drops in a periodic shear flow［J］.Journal of Non－Newtonian Fluid Mechanics，2004，118（1）：65.

［11］ Boek E S，Coveney P V，Lekkerkerker H N W.Computer simulation of rheological phenomena in dense colloidal suspensions with dissipative particle dynamics［J］.J Phys：Condens Matter，1996，8：9509.

［12］ Fan X J，Phan－Thien N，Chen S，et al.Simulating flow of DNA suspension using dissipative particle dynamics［J］.Physics of Fluids，2006，18（6）：063102.

［13］ Guo X D，Tan J P K，Zhang L J，et al.Phase behavior study of paclitaxel loaded amphiphilic copolymer in two solvents by dissipative particle dynamics simulations［J］.Chemical Physics Letters，2009，463：336.

［14］ Boek E S，Coveney P V，Lekkerkerker H N W，et al.Simulating the rheology of dense colloidal suspensions using dissipative particle dynamics［J］.Physical Review E，1997，55：3124.

［15］ Espa?ol P，Waren P B.Statistical mechanics of dissipative particle dynamics ［J］. Europhysics Letters，1995，30（4）：191.

［16］ Huilgol R R，Phan－Thien N.Fluid mechanics of viscoelasticity：general principles，constitutive modeling，analytical and numerical techniques［M］.Amsterdam：Elsevier，1997.

［17］ Fan X J，Phan－Thien N，Ng T Y，et al.Microchannel flow of a macromolecular suspension［J］.Physics of Fluids，2003，15（1）：11.

［18］ Rapaport D C.The art of molecular dynamics simulation［M］.New York：Cambridge University Press，2004.

［19］ Groot R D，Warren P B.Dissipative particle dynamics：bridging the gap between atomistic and mesoscopic simulation［J］.Journal of Chemical Physics，1997，107（11）：4423.

［20］ Evans D J， Morriss G P. Statistical mechanics of nonequilibrium liquids［M］.Canberra：ANU E Press，1990.

［21］ Chatterjee A. Modification to Lees－ Edwards periodic boundary condition for dissipative particle dynamics simulation with high dissipation rates［J］.Molecular Simulation，2007，33（15）：1233.