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

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

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

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

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

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

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

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

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

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

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

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

式中的求和符號對一定的截斷半徑rc內(nèi)所有的不包括粒子i的其他粒子進(jìn)行求和.如果粒子之間的距離rij＞rc，則粒子之間的相互作用為零.

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

式中：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時，wD＝0；vij為粒子之間相對速度，vij＝vi－vj.隨機(jī)力FR，ij可表達(dá)為

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

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

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

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

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

1.2 粒子間作用力的計算

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

由于粒子間相互作用力的對稱性，因此可以僅考慮相鄰的元胞中的一半數(shù)目，對于二維問題來說，是5個元胞（含中心元胞本身）；對于三維問題來說，是14個元胞（同樣也含中心元胞本身）.對于周期性邊界條件，相鄰的鏡像元胞也必須計算在內(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］來更新每個粒子在新時刻（即t＋Δt時刻）的位置、速度以及受力情況.如下所示：ri（t＋Δt）＝式中：為粒子i在t＋Δt時刻的預(yù)估速度；λc為經(jīng)驗系數(shù).如果粒子的受力情況與粒子間的相對速度無關(guān)，則λc＝0.5.Groot和Warren［19］發(fā)現(xiàn)，DPD方法中λc的最優(yōu)值為0.65.當(dāng)采用該最優(yōu)值同時取ρ＝3，σ＝3時，對平衡系統(tǒng)的模擬中計算的時間步長可以增加到Δt＝0.06而不失去對溫度的控制.

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

2.1 LE邊界條件

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

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

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

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

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

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

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

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

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

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

3 結(jié)果及討論

3.1 LE邊界條件下DPD簡單剪切流動模擬

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

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

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

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

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

Schmidt數(shù)Sc是動量擴(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ù)的一種有效途徑，從而使所模擬的流體更接近實際流體［12］.本文中模擬了σ＝14.14，即γ＝100情況下，時間步長為0.02、計算區(qū)域為30.00×10.00×22.25時系統(tǒng)的速度分布.如圖7所示，在較高施密特數(shù)條件下，剪切率為0.2時，系統(tǒng)內(nèi)的速度分布仍然是線性的，流動沒有邊界效應(yīng)，說明速度梯度在整個空間區(qū)域中是均勻的.

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

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

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

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

4 結(jié)論

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

［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)洪夫.耗散粒子動力學(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.