專利名稱：分子動力學模擬中壁面邊界的模擬方法
技術領域：
本發明涉及一種分子動力學模擬中的方法，尤其是模擬粒子在壁面邊界條件下的 運動軌跡的方法。
背景技術：
現在越來越多的科學家用分子動力學方法模擬平行平板中流體和納米通道內流 體流動的動力學特性。隨著研究的深入，科學家們開始對非平直納米通道進行研究和采用 分子動力學方法對此通道進行模擬，如=Xi-Jun Fan等人采用非平衡分子動力學方法模擬 了簡單流體在周期性噴管中的流動特性；為研究粗糙度對微通道內體流動及其邊界滑移性 質的影響，曹柄陽等人研究了氬氣在鋸齒狀鉬通道內的流動J. Castillo-Tejas等人分別 對線性聚合鏈和牛頓流體在4 1 4收縮-擴展通道中的流動過程進行了模擬。在使用 分子動力學方法模擬兩平板間和納米通道內流體的流動特性時，以前的文獻中對虛擬熱壁 面(Virtualthermal wall)情況下壁面邊界條件的模擬有兩種方法，一種是Allen，M. P. and D. J. Tildesley，在文獻Computer Simulation of Liquids.中提出的全反射法，全反射法 將壁面看成光滑鏡面，粒子在壁面上發生鏡面反射，采用全反射法模擬粒子在壁面邊界條 件下的運動軌跡參見圖1(a)，當粒子從A點運動到B’時，由于虛擬壁面的存在，粒子不會 在B’的位置出現，粒子會在壁面發生全反射后運動到B點，B與B’點關于壁面對稱，與壁 面平行坐標方向的速度矢量保持不變，與壁面垂直坐標方向的速度矢量大小保持不變，方 向反向；雖然全反射法可以真實地反應出粒子在發生發射后的位置，但是全反射法不能模 擬出粒子在定壁溫和壁面不光滑的條件下的運動軌跡，并且如果在非平直納米通道中，若 采用全反射法來模擬粒子的壁面邊界條件下的軌跡，則在計算速度矢量方面也存在一定的 X隹度；另一禾中是 Raraport，D. C.在文獻 The art of molecular dynamics simulation.中 提出了隨機反射法，隨機反射法與物理上的漫發射對應，采用隨機反射法模擬粒子在壁面 邊界條件下的運動軌跡參見圖1(b)，粒子反射到B’向壁面做垂線的垂點B，速度大小由壁 面溫度確定，方向隨機。雖然隨機反射法可以模擬出壁面不光滑以及定壁溫的物理條件，但 是隨機反射法由于粒子反射后的位置在壁面垂直方向的單一性，因此在模擬粒子運動軌跡 的過程中存在一定的不真實性，不能真實反映出粒子反射后的位置。因此，如何真實模擬粒子在壁面邊界條件下的運動狀態，成為了本領域內的研究 人員急需解決的技術問題。
發明內容
針對現有技術存在的上述不足，本發明提供一種在分子動力學模擬中，當粒子處 于壁面邊界條件下時，可綜合全反射法和隨機反射法優點，真實體現粒子在壁面邊界條件 下的運動狀態的壁面邊界模擬方法。為了實現上述發明目的，本發明采取以下技術方案分子動力學模擬中壁面邊界 的模擬方法，其特征在于，所述模擬方法包含以下步驟
A)定義分子動力學模擬中粒子的數據結構，粒子在三維模型中由六個信息來表示 粒子在分子動力學模擬通道中的運動狀態，所述六個信息分別是X方向坐標、y方向坐標、Z 方向坐標、χ方向速度、y方向速度和ζ方向速度；其中χ方向坐標為壁面橫坐標，y方向坐 標為壁面縱坐標，ζ方向為與XY平面垂直的方向；B)確定粒子在通道中的當前位置A(X(1，y0, z0, Utl，V0, W(1)，其中X(1表示粒子在χ方 向上的當前位置坐標，ytl表示粒子在y方向上的當前位置坐標，Ztl表示粒子在ζ方向上的 當前位置坐標，Utl表示粒子在χ方向上的當前位置速度，Vtl表示粒子在y方向上的當前位 置速度，W0表示粒子在ζ方向上的當前位置速度；C)根據粒子當前位置及粒子受到其它粒子的作用力信息，確定粒子的運動方向， 并根據粒子當前位置及受力信息判斷粒子下一步是否會運動到壁面邊界以外位置；如果粒 子沒有運動到壁面邊界以外位置，則粒子會沿著確定的運動方向運動，到達下一步指定位 置B’(x’，y’，Z’，U’，V’，w’)，其中χ’為粒子在χ方向上的指定位置坐標，y’表示粒子在 y方向上的指定位置坐標，ζ’表示粒子在ζ方向上的指定位置坐標，u’表示粒子在χ方向 上的指定位置速度，V’表示粒子在y方向上的指定位置速度，W'表示粒子在Z方向上的指 定位置速度；D)若粒子運動到壁面邊界以外位置，粒子將不能到達下一步指定的位置B’，粒子 根據壁面的溫度，經過壁面反射到達下一步反射位置B (X，y, z, U, ν, w)，當χ方向存在壁 面，則其中 y = y’，ζ = ζ’，χ = 2* 壁面 χ 方向的坐標-χ，，u=a* ^JkeTfm ,v=a* ^kBT/m，
Vf= a* Hm，粒子的速度方向隨機；當y方向存在壁面，則其中χ = χ’，ζ = y = 2* 壁面y方向的坐標_y’，VV> ^Tjm，w= a* ,]kBT/m，粒子的速度方向隨機； 當ζ方向存在壁面，則其中χ = x’，y = y’，z = 2*壁面ζ方向的坐標-ζ，，u=a* ^kBT/m， v=a*」kBT/m ,w= a* ^kBT/m，粒子的速度方向隨機。上述公式中，a為-1 1之間的隨機數(該隨機數的取得采用文獻The art of moleculardynamics simulation(D· C· Rapaport，Second Edition，Cambridge University Press, 2004)中的程序進行計算，在此不作詳細描述)，kB為波爾茲曼常數，T為當前溫度，m 為質量，χ表示粒子在χ方向上的反射位置坐標，y表示粒子在y方向上的反射位置坐標，ζ 表示粒子在ζ方向上的反射位置坐標，u為粒子在χ方向上的反射位置速度，ν為粒子y方 向上的反射位置速度，w為粒子在ζ方向上的反射位置速度。本發明分子動力學模擬中壁面邊界的模擬方法，其主要原理是在鏡面反射基礎上 考慮了壁面粗糙度對粒子速度的影響，因此，本發明方法又可以稱為半反射法。通常情況下，統計系統的宏觀量，可以通過粒子的運動狀態，即在本發明方法通過 粒子分別在X，y，Z方向上的坐標和粒子分別在X，y，Z方向的速度來反映粒子的運動狀態， 當粒子運動狀態發生變化時，χ, ι, ζ方向上的坐標和粒子分別在χ，y，ζ方向的速度都會發 生變化。當采用半反射方法對粒子處于壁面邊界條件下的運動軌跡進行模擬時，其既考慮 粒子的隨機性，即粒子在遇到壁面進行反射后，其反射位置速度(u，ν, w)會根據壁面的溫 度有關，而其反射位置x，y，ζ方向坐標(x，y，z)則與粒子未碰到壁面的預定的下一步指定 位置的χ，ι, ζ方向坐標(χ’，y’，ζ’ )有關。因此本方法在模擬粒子下一步位置時，不僅考 慮到壁面的溫度，也考慮到了粒子經過壁面發生反射后，其位置必定與粒子的未經過反射的位置有關系。本發明產生的有益效果(1)本發明方法不僅可以在分子動力學模擬中體現恒壁溫條件和壁面不光滑的物 理條件，而且也可以在分子動力學模擬中真實反映粒子反射后的位置。其解決了隨機反射 法中不能真實反應粒子反射后位置的問題。(2)本發明方法還可以在不規則納米通道中，簡化速度矢量計算方法，從而便于計 算。解決了全反射邊界法不能模擬恒壁溫條件以及壁面不光滑條件下粒子的真實運動狀態 的問題和在不規則納米通道中，全反射法不能簡化速度矢量計算方法的問題。(3)本發明方法可以用于模擬粒子在各種不同通道中運動到虛擬壁面邊界外之后 的運動情況，從而為研究粒子在不同通道中運動狀態，提供了技術支撐。
圖1是分子動力學模擬中的粒子運動軌跡示意圖；其中(a)為采用全反射法模擬 粒子在壁面邊界條件下的運動軌跡示意圖；其中(b)為采用隨機反射法模擬粒子在壁面邊 界條件下的運動軌跡示意圖；圖2是分子動力學模擬中的粒子在兩種邊界條件下反射示意圖其中(a)半反射 粒子軌跡運動示意圖；(b)隨機反射粒子軌跡運動示意圖；圖3是分子動力學模擬中的變截面通道示意圖；其中(a)為三維變截面通道示意 圖；(b)為三維變截面通道正視示意圖；圖4是分子動力學模擬中粒子數密度分布隨到壁面距離的比較圖；圖5是本發明方法模擬變截面納米通道中勢能隨時間變化曲線圖。
具體實施例方式下面結合附圖和具體實施方式
對本發明作進一步詳細地說明。本文中提出一種綜合了全反射法和隨機反射法優點的新的壁面邊界條件處理方 法——半反射法。當粒子A點運動到B’時，如圖1(a)，由于壁面的存在，粒子經過壁面反射 后到達B點，坐標位置與全反射法粒子反射后相同，但是速度大小由壁面溫度確定，方向隨 機，即速度矢量的確定方法與隨機反射法相同。在模擬非平直微納米通道時，比如非45° 傾斜通道，如果壁面采用全反射邊界條件處理法，粒子反射后，速度矢量的確定有一定的難 度。隨機反射法相對于全反射法隱含了兩個條件(1)壁面恒溫；(2)壁面不光滑，與真實 情況更加接近。但是隨機反射法在模擬納米通道時有一定的不真實性。如圖2所示，χ方 向為流動方向，當粒子運動到區域A時，如果采用隨機反射法，見圖2 (b)，無論粒子B處在A 中的任意位置經過坐標調節后都出現在B’，不能真實反應出粒子運動的坐標位置。但如果 采用半反射邊界條件模擬方法，如圖2(a)，粒子B出現在區域A中的任意位置，采用半反射 法經過坐標調節后可以出現在B對應任意B’點，從而與粒子的真實運動位置更為接近。實施例1本實施例通過在分子動力學模擬過程中，分別采用隨機反射法、全反射法和半反 射三種壁面邊界條件模擬方法來模擬平行平板中的氬流體流動過程，從而進一步說明本發 明方法與現有的隨機反射法和全反射法的不同之處。具體實施如下
采用氬原子模擬NVT系中，用速度修正法控制系統溫度，原子間相互作用采 用12-6LJ勢能模型，運動方程采用Verlet數值積分法，時間步長為lfs，模擬時間為 lOOOOOfs，模擬尺寸為6nmX6nmXl. 5nm，粒子數為648，密度為lg/cm3.在本實施例中，將 坐標原點定義為0點，χ和y方向采用周期性邊界條件，當ζ方向(Z方向是指三維通道中 與xy平面垂直的方向，粒子只可能超出ζ方向的壁面邊界，因為χ和y方向采用周期性邊 界條件，所以在這個實施例中都是無限延伸的，沒有壁面)為0時，xy平面為下壁面，當ζ 為1.5nm時，xy平面為上壁面，下面分別采用三種不同的壁面邊界條件模擬方法來獲取粒 子的運動狀態隨機反射法χ= χ，；y = y，；z = 0 ；M= a*^jk13T/m ；V= a*^kBT/m ；w= a*^jkeT/m全反射法x= x，;y = y，;z = 2* 壁面 ζ 方向的坐標 _z，;u = u，;v = v，;w = i’半反射法x= X，；y = y，；z = 2* 壁面 Z 方向的坐標-Z，-,U= a*^keTIm ； ν= a*yjkBT/m ；w= a*yJkBT/m上述公式中，a為-1 1之間的隨機數，kB為波爾茲曼常數，T為當前溫度，m為 質量，X表示粒子在X方向坐標上的反射位置坐標，y表示粒子在y方向坐標上的反射位置 坐標，ζ表示粒子在ζ方向坐標上的反射位置坐標，u為粒子在χ方向坐標上的反射位置速 度，ν為粒子在y方向坐標上的反射位置速度，w為粒子在ζ方向坐標上的反射位置速度。 本文中的“*”為乘號，“_”為減號。從采用隨機反射法、全反射法和半反射法的模擬在特定條件下平行平板中的氬 流體流動中碰壁的情況來看，當采用隨機反射法模擬時，參見圖1(b)粒子A(X(I，y0, z0, u0, v0, w0)在運動過程中超出模擬區域邊界時，粒子不會沿著既定線路到達下一步指定位 置B’(χ’，y’，ζ’，U’，V’，W’)，而會沿著粒子反射到B’向壁面做垂線的垂點B，χ = χ’ ；y =y’ ；ζ =壁面ζ方向的坐標；速度大小由壁面溫度確定，《=口*彳‘7>2 ,ν= a*^kJjm、 W= a*^bTIm ；方向隨機，；當采全反射法模擬式，參見圖1(a)，粒子A (X(1，％，、，％，％，％) 在運動過程中超出模擬區域邊界時，粒子不會沿著既定線路到達下一步指定位置B’(χ’， y，，ζ，，U，，V，，W，)，粒子會在壁面發生全反射后運動到B (x, y, ζ, u, ν, w)點，B (χ, y, ζ, u, v,w)與B’(X’，y’，Z’，U’，V’，w’)點關于壁面對稱，與壁面平行坐標方向的速度矢量保持 不變，與壁面垂直坐標方向的速度大小保持不變，速度方向反向；當采用半反射法模擬時， 粒子A(X(I，y0, z0, u0, v0, w0)在運動過程中超出模擬區域邊界時，粒子不會沿著既定線路到 達下一步指定位置B’(χ’，y’，ζ’，U’，V’，W’)，而會在虛擬的壁面邊界上發生反射，到達反 射位置B (x, y, z,u, ν, w)，其反射后位置B既與粒子未發生反射時的指定位置B’(χ’，y’， z’，u’，v’，W’)有關，也與粒子發生反射時的壁面的溫度有關，即χ = χ’ ；y = y’ ；z = 2* 壁面ζ方向的坐標-ζ， ’u= a*Hm -,ν= a*^keTlm ；w= a*^kBT/m，粒子反射后的速度 方向隨機，從而得到粒子反射后的位置。由于粒子反射后的速度與三個方向的速度分量u， ν和w有關，因此，粒子反射后的速度大小為λΑ/2 +V2 + V ,而速度的方向為三個速度分量的 速度方向的合成，而速度分量是有正有負的，所以合成的速度方向也是各個角度都有的，從 而粒子反射后的速度方向是隨機的。本說明書中，壁面χ方向的坐標、壁面y方向的坐標和 壁面ζ方向的坐標是指壁面在三維笛卡爾坐標系的坐標，由于在實際應用過程中三維坐標系的原點可以為0點，也可以不為0點，因此，壁面本身在三維坐標系中的坐標與三維坐標 系原點的設置有關。2)模擬結束后，統計了三種壁面邊界條件模擬方法情況下粒子數密度分布的情 況，如附圖4所示。比較結果說明采用半反射邊界條件模擬法和采用全反射邊界條件模擬法模擬得 到的結果100%吻合，具有很強的實用性。從而表明本發明方法在模擬粒子在壁面邊界條件 下的運動狀態時，不僅可以體現恒壁溫條件和壁面不光滑的物理條件，而且也可以在分子 動力學模擬中真實反映粒子反射后的位置。半反射邊界法與隨機反射邊界法和全反射邊界法異同見表1和表2 表1半反射邊界法與隨機反射邊界法和全反射邊界法特性比較表
權利要求
分子動力學模擬中壁面邊界的模擬方法，其特征在于，所述模擬方法包含以下步驟A)定義分子動力學模擬中粒子的數據結構，粒子在三維模型中由六個信息來表示粒子在分子動力學模擬通道中的運動狀態，所述六個信息分別是x方向坐標、y方向坐標、z方向坐標、x方向速度、y方向速度和z方向速度；其中x方向坐標為壁面橫坐標，y方向坐標為壁面縱坐標，z方向為與xy平面垂直的方向；B)確定粒子在通道中的當前位置A(x0，y0，z0，u0，v0，w0)，其中x0表示粒子在x方向上的當前位置坐標，y0表示粒子在y方向上的當前位置坐標，z0表示粒子在z方向上的當前位置坐標，u0表示粒子在x方向上的當前位置速度，v0表示粒子在y方向上的當前位置速度，w0表示粒子在z方向上的當前位置速度；C)根據粒子當前位置及粒子受到其它粒子的作用力信息，確定粒子的運動方向，并根據粒子當前位置及受力信息判斷粒子下一步是否會運動到壁面邊界以外位置；如果粒子沒有運動到壁面邊界以外位置，則粒子會沿著確定的運動方向運動，到達下一步指定位置B’(x’，y’，z’，u’，v’，w’)，其中x’為粒子在x方向上的指定位置坐標，y’表示粒子在y方向上的指定位置坐標，z’表示粒子在z方向上的指定位置坐標，u’表示粒子在x方向上的指定位置速度，v’表示粒子在y方向上的指定位置速度，w’表示粒子在z方向上的指定位置速度；D)若粒子運動到壁面邊界以外位置，粒子將不能到達下一步指定的位置B’，粒子根據壁面的溫度，經過壁面反射到達下一步反射位置B(x，y，z，u，v，w)，當x方向存在壁面，則其中y＝y’，z＝z’，x＝2*壁面x方向的坐標 x’，粒子的速度方向隨機；當y方向存在壁面，則其中x＝x’，z＝z’，y＝2*壁面y方向的坐標 y’，粒子的速度方向隨機；當z方向存在壁面，則其中x＝x’，y＝y’，z＝2*壁面z方向的坐標 z’，粒子的速度方向隨機；上述公式中，a為 1～1之間的隨機數，kB為波爾茲曼常數，T為當前溫度，m為質量，x表示粒子在x方向上的反射位置坐標，y表示粒子在y方向上的反射位置坐標，z表示粒子在z方向上的反射位置坐標，u為粒子在x方向上的反射位置速度，v為粒子在y方向上的反射位置速度，w為粒子在z方向上的反射位置速度。FSA00000292915600011.tif,FSA00000292915600012.tif,FSA00000292915600013.tif,FSA00000292915600014.tif,FSA00000292915600015.tif,FSA00000292915600016.tif,FSA00000292915600017.tif,FSA00000292915600018.tif,FSA00000292915600019.tif
全文摘要
本發明公開提出一種分子動力學模擬中壁面邊界的模擬方法，該方法綜合了全反射法和隨機反射法優點，能夠真實反應出粒子碰壁后的運動狀態，也稱為半反射法，其實現步驟如下1)定義分子動力學模擬中粒子的數據結構；2)確定粒子在通道中的當前位置A(x0，y0，z0，u0，v0，w0)；3)根據粒子當前位置及粒子受到其它粒子的作用力信息，確定粒子的運動方向，并根據粒子當前位置及受力信息判斷粒子下一步是否會運動到壁面邊界以外位置；4)若粒子運動到壁面邊界以外位置，經過壁面反射到達下一步反射位置B(x，y，z，u，v，w)，當x方向存在壁面，則其中y＝y’，z＝z’，x＝2*壁面x方向的坐標-x’，粒子的速度方向隨機。本發明在模擬中可以體現出壁面恒溫和壁面不光滑的物理條件；可以更加真實地反應出粒子反射后的位置；也可應用在不規則的納米通道中，簡化速度矢量的計算方法。
文檔編號G06F17/50GK101944151SQ201010299639
公開日2011年1月12日 申請日期2010年9月30日 優先權日2010年9月30日
發明者丁玉棟, 馮潔, 葉丁丁, 吳睿, 廖強, 朱恂, 李俊, 王宏, 王永忠 申請人:重慶大學