study on microscope characteristics of interface by molecular dynamics simulation a thesis submitted to chongqing university in partial fulfillment of the requirement for the degree of doctor of engineering by cai zhiyong supervisor: prof. zeng danling major: power engineering 表 2.13 列出了甲烷在溫度 t*1.015 條件下，模擬的初始條件,表中模擬盒的長度 l*由公 式(2.16)計算得出： 表 2.12 氬氣過熱狀態(tài)的模擬體系初始條件(n=1372) table 2.12 initial conditions of molecular dynamics simulation for super-heated vapor of argon v*/v0* * l* 1 0.0611 28.21 2 0.0306 33.54 4 0.0153 44.78 8 0.0076 56.42 2 實際流體 p-v-t 特性的分子動力學模擬研究 33 表 2.13 甲烷氣體過熱狀態(tài)模擬體系初始條件(n=1372) table 2.13 initial conditions of molecular dynamics simulation for super-heated vapor of methane v*/v0* * l* 1 0.03408 34.27 2 0.01704 43.18 4 0.00852 54.40 8 0.00426 68.54 2.4.2 計算結果及分析 計算結果 首先， 將氬氣和甲烷看作理想氣體， 利用無量綱理想氣體狀態(tài)方程式(2.24)式， 可以計算得到在相同溫度、不同體積膨脹比條件下，氬氣和甲烷的過熱蒸氣壓力 值;然后， 將氬氣和甲烷看作實際氣體， 利用無量綱二常數半經驗實際氣體狀態(tài) p-r 方程式(2.29)式，計算可得到在相同溫度、不同體積膨脹比條件下，氬氣和甲烷的 過熱蒸氣壓力值; 最后，采用分子動力學模擬方法，利用式 (2.18)、(2.19)可以得 到不同初始條件氬氣和甲烷模擬體系的過熱蒸氣壓力。表 2.14、表 2.15 分別列出 了氬氣在溫度 t*1.088 和甲烷在溫度 t*1.015 的條件下，采用上述三種方法計 算的體積 v*v0*、2v0*、4v0*、8v0*時的過熱蒸氣壓力。 表 2.14 氬氣不同體積膨脹比的過熱蒸氣壓力計算值對比(t*=1.088) table 2.14. the comparison of pressure of the super-heated vapor at different volume ratio of expansion for ar(t*=1.088) v*/v0* 理想氣體 p-r md1 md2 1 0.0666 0.0485 0.0490 0.0487 2 0.0333 0.0285 0.0291 0.0288 4 0.0167 0.0154 0.0154 0.0159 8 0.0083 0.0080 0.0082 0.0083 表 2.15 甲烷不同體積膨脹比的過熱蒸氣壓力計算值對比(t*=1.015) table 2.15. the comparison of pressures of the super-heated vapor at different volume ratio of expansion for ch4(t*=1.015) v*/v0* 理想氣體 p-r md1 md2 1 0.0346 0.0280 0.0284 0.0276 2 0.0173 0.0156 0.0161 0.0158 4 0.0087 0.0082 0.0083 0.0086 8 0.0043 0.0042 0.0043 0.0044 重慶大學博士學位論文 34 圖 2.7 氬氣在不同體積膨脹比的過熱蒸氣壓力 fig.2.7 pressures of the super-heated vapor at different volume ratio of expansion for ar 表中理想氣體項是將氬氣和甲烷看作理想氣體計算得到的壓力值；表中 p-r項則是將氬氣和甲烷看作實際氣體，計算得到的壓力值; md1 表示采用分 子動力學模擬方法具體第一種算法計算的壓力，md2 表示采用分子動力學模擬方 法具體第二種算法計算的壓力。圖 2.7、圖 2.8 是相應的示意圖。 圖 2.8 甲烷在不同體積膨脹比的過熱蒸氣壓力 fig.2.8 pressures of the super-heated vapor at different volume ratio of expansion forch4 0123456789 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 p * v */v 0 * ideal gas p-r md1 md2 t *=1.088 0123456789 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 p * v */v 0 * ideal gas p-r md1 md2 t*=1.015 2 實際流體 p-v-t 特性的分子動力學模擬研究 35 結果分析 表 2.16、表 2.17 分別列出了采用理想氣體模型和分子動力學模擬方法(md1、 md2)計算的氬氣和甲烷的過熱蒸氣壓力與氬氣和甲烷的過熱蒸氣壓力的實驗值 相比的相對誤差。需要說明的是：由于缺乏氬氣和甲烷的過熱蒸氣壓力實驗值數 據，在下面的比較中，我們將采用實際氣體 p-r 方程計算的氬氣和甲烷的過熱蒸 氣壓力值作為參照值進行對比。從上節(jié)中的計算可以看出，在飽和狀態(tài)下，p-r 方 程計算的氬氣和甲烷的飽和壓力值與其參照值相差不到 3%；在此基礎上，體積進 一步等溫膨脹到過熱狀態(tài)，預計采用 p-r 方程計算的氬氣和甲烷的過熱蒸氣壓力 值與之相比，誤差會進一步減小，因此，用 p-r 方程計算的氬氣和甲烷的過熱蒸 氣壓力值替代過熱蒸氣壓力實驗值，從技術上講是合理的、可行的。 表 2.16 氬氣在不同膨脹比下的壓力計算值與實驗值的相對誤差(%) tab.2.16 the relative error (%)of pressure of the super-heated vapor at different volume ratio of expansion for ar(t*=1.088) v*/v0* 1 2 4 8 理想氣體 37.3 16.8 8.4 3.75 md1 0.71 3.21 1.22 2.38 md2 0.41 1.05 2.60 3.75 表 2.17 甲烷在不同膨脹比下的壓力計算值與實驗值的相對誤差(%) tab. 2.17 the relative error (%)of pressure of the super-heated vapor at different volume ratio of expansion for ch4(t*=1.015) v*/v0* 1 2 4 8 理想氣體 23.57 12.34 6.10 2.38 md1 1.43 3.21 1.22 2.38 md2 1.43 1.28 4.88 4.76 從表 2.16、表 2.17 可以看出： 1)采用分子動力學模擬方法計算的氬氣和甲烷的過熱蒸氣壓力與用 p-r 方程 所得數據相比，結果接近，其相對誤差均不超過 5%。 2)采用理想氣體模型計算的氬氣和甲烷的過熱蒸氣壓力與實驗值相比， 相對誤 差最大，但隨著膨脹比的逐步增大，其與實驗值相比的相對誤差逐步減小，當膨 脹比大于 8 后，無論對氬氣或是對甲烷，其相對誤差均不超過 5%，可以采用理想 氣體模型直接計算氬氣或對甲烷的 p-v-t 特性。 2.5 實際氣體特性的分數布朗運動解釋 從上面兩節(jié)的模擬結果可以總結出以下三個結論： 重慶大學博士學位論文 36 對 l-j 流體，在計算實際氣體的 p-v-t 特性上，md 方法所得數據具有較 高精度，其數值級與用和 p-r 狀態(tài)方程式所得數據相當；因此，md 方法可近似用 于 l-j 流體飽和狀態(tài)的熱力學特性計算。 理想氣體模型，由于它是無相變的理想流體模型，不能用于計算實際氣體 (即使是 l-j 流體)飽和狀態(tài)或近飽和狀態(tài)的熱力學特性計算； 特別是飽和溫度越高， 采用理想氣體模型計算的誤差越大。 當實際氣體遠離飽和狀態(tài)時(本文中實際氣體比容大于 8 倍同溫度下的飽 和比容)，可以近似采用理想氣體模型計算實際氣體的 p-v-t 特性。 運用熱力學知識，可以對這種現(xiàn)象進行合理的解釋：當氣體處于飽和狀態(tài)時， 隨著溫度升高，飽和壓力也升高，壓力越高，實際氣體與理想氣體偏差越大;同樣, 當氣體在飽和狀態(tài)下等溫膨脹時，隨著膨脹比增加,氣體壓力降低，壓力越低，實 際氣體與理想氣體偏差越小。 運用分子動力學模擬方法我們還可以從定量上解釋這種現(xiàn)象，以重慶大學曾 丹苓教授為首的課題組對造成這種現(xiàn)象的原因進行了深入的研究，曾丹苓教授在 mendelbrot 將布朗運動的有關理論擴大到分數布朗運動119 120的基礎上進一步提 出用分數布朗函數作為描寫實際流體分子無規(guī)則運動的概率密度函數，從而從根 本上解釋了實際氣體與理想氣體產生誤差的原因，將分形理論運用于分子隨機運 動的基本論點可以歸納如下： 首先，標準布朗運動是一隨機過程，粒子的位移 x 是一隨機變量，滿足： 函數 x(t)連續(xù)，且以概率為 l 滿足 x(0)=0; 對任意 t0 及 h0， 位移增量滿足均值為零， 方差為 h 的正態(tài)分布， 即 x(t+h) x(t)n(0，h)，其相應的分布函數為： du h u h xtxhtxp x 2 exp 2 1 2 (2.33) 概率密度函數為： e h x h tx 2 2 2 1 , (2.34) 若 0t1t2tm，位移增量 x(t2)-x(t1)，x(tm)-x(tm-1)彼此獨立。 由以上定義可見，位移 x 的數學期望 e(x)為零，而方差為： dxxxxe b 22 2 (2.35) 其次，mendelbrot 在 1968 年，將布朗運動的有關理論擴大到分數布朗運動， 并給出以下定義 2 實際流體 p-v-t 特性的分子動力學模擬研究 37 分數布朗運動是隨機過程，粒子位移 x 滿足： 函數 x(t)連續(xù)，且以概率為 l 滿足 x(0)=0； 對任意 t0 及 h0，位移增量滿足均值為零，方差為 h2h的分布，即 x(t+h) x(t)n(0，h2h)，其分布函數為： du h u h xtxhtxp x h h 2 2 2 2 exp 2 1 (2.36) 概率密度函數為 e h h x h h tx 2 2 2 2 2 1 , (2.37) 若 0tlt2tm，位移增量 x(t2)-x(t1)，x(tm)-x(tm-1)相互關聯(lián)。 這里，h(0，1)稱為 hurst 指數。當 h1/2 時，以上的隨機過程即定義為分數布 朗運動，而當 h1/2 時，則還原為標準布朗運動。可以證明，當 0h1 時，以上 的定義仍成立。 從分形理論可知，圖像的分數維數用下式求得121 d=3-h (2.38) 式中，h 為 hurst 指數，3 為空間維數。其中，維數 d 不一定為整數，常包含整數 和分數兩個部分，整數部分表示運動的拓撲維數，而分數部分則用來度量無規(guī)運 動的復雜程度及其充滿空間的能力。 嚴格地說，標準布朗函數僅適用于描述處于平衡態(tài)的理想氣體的分子運動。 對于分子運動具有分形特征的實際流體，系統(tǒng)中分子、原子等微粒的無規(guī)則運動 可以用分數布朗函數作為概率密度函數來描寫。 同時，構造出以下形式的微粒位移 x 的分布函數： du kh u kh xtxhtxp x h h 2 2 2 2 exp 2 1 (2.39) 相應的概率密度函數為： h h kh x kh txf 2 2 2 2 exp 2 1 ),( (2.40) 其中，位移 x 的平均仍為零，方差為 h khtxhtxe 22 2 (2.41) 式中，h 為時間增量；k、h 分別為待定的系數和指數，可根據粒子運動的圖 像確定，粒子運動的圖像一方面可從實驗觀測，另一方面可通過計算機模擬得出。 重慶大學博士學位論文 38 將上式對上式兩邊取對數得： khhlnln2ln 2 (2.42) 這樣，待定系數 k 和指數 h 可由hlnln 2 的關系求得，在對數坐標hlnln 2 上 曲線的斜率即是 2h，而在縱軸上的截距即是hln。 在分子動力學模擬中，位移 x 的方差可直接根據數值模擬的結果求得，曾丹 苓教授的博士生劉娟芳69運用分子動力學模擬得到的位移的數據，采用最小二乘 法進行擬合，得到了實際流體 ar 不同溫度下的 hurst 指數和分維數 d。 本文采用 同樣的方法計算了氬氣在溫度 t*1.088 和甲烷在溫度 t*1.015 的條件下，體積 v*v0*、2v0*、4v0*、8v0*時的 hurst 指數和分維數 d，表 2.18 列出了不同膨脹 比下氬氣的 hurst 指數 h 和分數維數 d 表 2.18 氬在不同體積膨脹比的 hurst 指數 h 及分維數 d(t*=1.088) table 2.18 the hurst exponential h and fractional dimensionality d at different volume ratios of expansion for ar(t*=1.088) v*/v0* 1 2 4 8 h 0.5155 0.5126 0.5098 0.5073 d 2.4845 2.4874 2.4902 2.4927 表 2.19 列出了不同膨脹比下甲烷的 hurst 指數 h 和分數維數 d： 表2.19 甲烷在不同體積膨脹比的hurst指數h及分維數d(t*=1.015) table 2.19 the hurst exponential h and fractional dimensionality d at different volume ratios of expansion for ch4(t*=1.015) v*/v0* 1 2 4 8 h 0.5125 0.5106 0.5098 0.5068 d 2.4875 2.4894 2.4902 2.4932 從表 2.18、2.19 可以看出，無論是氬氣或是甲烷，其分數維數 d 及 hurst 指 數 h 與氣體等溫膨脹的程度有關。等溫膨脹體積比越大，氬氣和甲烷的 hurst 指 數越加接近理想氣體的 hurst 指數 0.5，因此，采用分數布朗運動我們可以描述實 際氣體偏離理想氣體的程度。 2.6 小結 在本章中, 主要作了三個工作: 2 實際流體 p-v-t 特性的分子動力學模擬研究 39 通過對單相 l-j 流體進行分子動力學模擬和實際氣體的 p-v-t 特性計算, 我 們得出結論: 用 md 方法可得到一定精度的熱力學數據其結果與用 p-r 實際氣體 狀態(tài)方程式計算的 l-j 流體的 p-v-t 特性比較甚為吻合，因此，md 方法用于 l-j 流體飽和狀態(tài)和過熱蒸氣狀態(tài)的熱力學特性計算是近似可行的。 本章除用簡單的 l-j 流體進行了飽和參數的計算外， 還將模擬計算擴大到較 為復雜，而在工程上用得較多的強極性分子水的飽和參數的計算，結果與實驗值 吻合極好，進一步證明了所用方法的可行性和其工程意義。特別是對一些缺乏飽 和參數實驗值的物質，模擬計算不失為一種有效的方法以補充其數據的不足，滿 足工程的需要。 運用將用分數布朗函數作為描寫實際流體分子無規(guī)則運動的概率密度函數 的理論，利用分子動力學模擬過程中產生的大量的分子運動的微觀信息，計算得 出了氬氣和甲烷從飽和狀態(tài)等溫膨脹到不同的過熱蒸氣狀態(tài)下的分數維數 d 及 hurst 指數 h，證明等溫膨脹體積比越大，氬氣和甲烷的 hurst 指數越加接近理想 氣體的 hurst 指數 0.5，可以用分數布朗運動我們描述實際氣體偏離理想氣體的程 度。 從本章可初步看出分子動力學模擬在研究物質物性方面的重要作用，以上的 分析可認為僅反映了冰山的一角， md 模擬方法還可擴大到更多的方面， 以下各章 將著重應用 md 于界面現(xiàn)象的分析，得到更多具有理論和工程意義的結果。 重慶大學博士學位論文 40 3 氣-液界面微觀特性的分子動力學模擬研究 41 3 氣-液界面微觀特性的分子動力學模擬研究 流體的分子動力學模擬近年來得到很大的發(fā)展，已證明是一個研究流體物性 的重要工具。運用分子動力學模擬方法研究氣-液界面現(xiàn)象，是近二十年來，隨著 gibbs 系綜的應用，極大地推動了相平衡的分子模擬技術發(fā)展的最成功的事例 49,50。 界面是不同相間的邊界面。由于界面非常薄，在許多場合下僅有幾個分子厚， 因而常被忽略而將其想像為一個沒有厚度的幾何面8,9，但正是這一薄薄的界面的 存在，對界面附近的流動和輸運、化學過程產生了重要的影響，例如蒸餾過程中 的輸運速度、液體的萃取、表面活化等工藝過程中其影響都是顯著的，特別是對 于當前迅速發(fā)展的微尺度物理過程，其影響更不容忽略。 本章的目的是從微尺度上借助分子動力學方法， 以最常見的氣-液界面層為例， 系統(tǒng)、全面地研究界面層的物理特性，包括其中分子的密度分布、溫度分布、表 面層厚度、表面張力及其影響因素。在模擬計算中，作者提出了一種新的變截斷 半徑算法以提高計算精度、改進計算過程，取得了較好效果。 采用分子動力學模擬方法對氣-液共存相進行模擬最直接的辦法是將模擬盒分 成明確的區(qū)域。例如，將氣-液相平衡放在包含有氣相區(qū)、液相區(qū)、氣-液界面區(qū)等 三個明確的區(qū)域的模擬盒中進行，模擬從將預先已平衡的液體塊放在模擬盒的中 央開始，如果液體的密度接近共存相的密度，液體將蒸發(fā)進入剩余空間形成穩(wěn)定 的兩相系統(tǒng)。 在本章里，將運用前面介紹的界面層模型和分子動力學模擬方法對液體和自 身的蒸汽形成氣-液界面的界面特性、界面張力、界面層中的密度、溫度等熱力學 參量的變化規(guī)律進行研究。 3.1 氣-液界面特性的分子動力學模擬研究方法 本文進行的氣-液界面特性研究，仍是以 l-j 流體為基礎，對 l-j 流體，分子 i 和 j 間的相互作用勢函數表示為1： 612 4)( ijij ij rr ru (3.1) 式中，為能量參數，為尺度參數，都是與具體流體性質相關的常數(氬和甲烷 的物性參數見表 2.1)， ij r為分子i和j之間的距離。 運用 lennard-jones 勢函數進行分子動力學模擬中，采用直角坐標系 0 xyz， 重慶大學博士學位論文 42 在初始時刻流體處于飽和液相，分子以面心立方晶格(fcc)均勻分布在邊長為l的 正方體模擬盒中，其質心在原點(0，0，0)處。模擬過程中，系統(tǒng)采用正則系綜， 即分子數n、容積v和溫度t保持不變，并不斷調整體系質心使之處于原點。模 擬盒在 x、y 方向的長度 lx=ly=l，為了保證足夠的氣相空間，取在z方向的長度 lz=3l(見圖 3.1)。本模擬在三個方向上都采用周期性邊界條件。 在計算機分子動力學模擬過程中，為確保模擬體系各物理量之間的數量級較 為接近， 減小模擬過程中出現(xiàn)的誤差， 利用 l-j12-6 勢函數中體系中的參數和， 以及分子質量m對模擬流體各物理量按照第二章表 2.1 進行無量綱化處理： 3.1.1 計算公式 在進行分子動力學模擬時， 分子的運動規(guī)律由聯(lián)立求解各分子的運動方程 2 2 dt rd mrff ij i ij iji ni 2, 1 (3.2) 得出，對 l-j 流體, 取截斷半徑為r時,分子間相互作用力 ij rf為1： ij ij ij dr rdu rf rr rr rrr ij ij ijijij , ,0 2 24 612 (3.3) 本文模擬的最終結果是計算氣-液相平衡中均勻相的飽和氣體、飽和液體的密 度以及非均勻相的表面張力， 模擬中將模擬盒沿z方向平行于x-y平面平均分為 s n 圖 3.1 模擬盒尺寸圖 fig.3.1 sketch of simulation geometry ly lz vapor vapor liquid lx 3 氣-液界面微觀特性的分子動力學模擬研究 43 個切片進行數理統(tǒng)計，每個切片的分子數密度： sl k v n k (3.4) 任一層k的溫度按下式計算： kn i i k k v n kt 1 2 3 1 (3.5) 界面層中非均勻相的法向應力、切向應力及表面張力122： ji ij ij ijk sl bn ru r z v tkkkp , 2 1 (3.6) ji ij ij ijij k sl bt ru r yx v tkkkp , 22 2 1 (3.7) ji ij ij ijijij k ru r zyx lxly k , 222 2 2 1 (3.8) s n k k 1 2 1 (3.9) 式(3.4)(3.7)中： k為第k切片的分子數密度， kt為第k切片的溫度， kpn、 kpt 為第k切片的法向應力和切向應力， k n為第k切片內的分子數，vi為分子i的速度， 式(3.8)中 k為第k切片的局域表面張力， ij x， ij y， ij z及 ij r分別為分子i和j之間在 x、y、z 方向及空間r的距離，vsl=lxlylz/ns為一個切片的體積，u為勢函數對 r 的導數，由式(3.3)計算得出， 為正則系綜統(tǒng)計平均值。在求和之中，如果i和j 兩個分子都在第k切片內，則將其全部計入；如只有其中一個分子位于第k切片內， 則計入其貢獻的一半。 3.1.2 氣-液界面特性分子動力學模擬中的技術問題 運用前面的分子動力學模擬方法，我們能夠對氣液兩相系統(tǒng)進行分子動力學 模擬，由于氣液兩相系統(tǒng)的特殊性，為了通過分子動力學模擬方法準確地反映氣- 液界面的宏觀和微觀特性，首先對影響分子動力學模擬結果的主要因素作重點的討 論。 模擬分子數和模擬時間問題 根據分子運動各態(tài)遍歷原理，模擬系統(tǒng)時空和分子數的不足可以通過延長系統(tǒng) 模擬時間來彌補，分子數不足帶來的尺度效應則可以通過選取合適的邊界條件彌 補。因此，從理論上講，模擬系統(tǒng)對分子數的多少并無嚴格要求，但實際上，模擬 重慶大學博士學位論文 44 的分子數過少，分子參數的波動帶來的影響就會加大，同時，分子數的多少決定了 模擬盒的尺寸，限制了勢能截斷半徑的選取，小的勢能截斷半徑選取將帶來較大的 計算誤差；反之，模擬的分子數也并非越多越好，模擬系統(tǒng)的分子數越多，達到平 衡所需的分子動力學步越少，但總的模擬時間不會減少，甚至會成幾何級數增加； m.mecke74等采用分子數 n=1372， 時間步長t*=0.002 對氣-液界面特性進行了分子 動力學模擬，證明均勻相、非均勻相達到平衡所需的時間是不同的，在均勻相要達 到平衡只需 1000020000 分子動力學步， 而在非均勻相達到平衡所需的時間則要長 得多。scott11等人對模擬分子數和模擬時間問題進行了系統(tǒng)的研究，表 3.1 是 scott 等人推薦的模擬分子數和模擬時間選取值。 表 3.1 模擬體系達到平衡時所需的時間 tab.3.1 equilibration time of simulation system 溫度 t* 系統(tǒng)的分子數 n 達到平衡的時間 0.8 7200 6000 24000 6000 69350 6000 0.9 14400 11000 40000 10000 1.0 14400 16500 48000 14400 注：達到平衡的時間單位為分子動力學步； 分子動力學步時間步長為 t* = 0.006。 從表 3.1 提供的研究結果可以看出， 系統(tǒng)達到平衡所需的分子動力學步不但受 到模擬分子數的影響，而且也受到系統(tǒng)溫度的影響，在分子數相同的條件下，系 統(tǒng)的溫度越高，系統(tǒng)達到平衡所需的分子動力學步數越多。 如何判斷氣-液系統(tǒng)是否達到平衡態(tài)？根據經典熱力學理論， 氣-液系統(tǒng)達到相 平衡，必須滿足三個平衡條件：即，1、熱平衡條件，系統(tǒng)內部的溫度均勻；2、 力平衡條件，系統(tǒng)內部的壓強均勻；3、相平衡條件，系統(tǒng)內各相的化學勢相等。 為了判斷分子動力學模擬過程中，系統(tǒng)何時達到相平衡條件，需要不斷地對系統(tǒng) 的能量、氣-液界面的密度分布、氣相和液相的分子數以及壓力分布等進行統(tǒng)計。 當分子動力學模擬系統(tǒng)達到平衡時，氣相區(qū)和液相區(qū)的分子數應基本保持不變， 即系統(tǒng)的密度分布和壓力分布基本上不再隨時間的變化而變化。 3 氣-液界面微觀特性的分子動力學模擬研究 45 本文對模擬分子數 n 對模擬時間的影響進行了專題研究，下表列出了勢能截 斷半徑r *=5.0，時間步長t*=0.001 下，不同分子數 n 運行 10000 分子動力學步 計算機所耗的計算時間(計算機為方正，cpu：piv；內存 256m)和達到氣-液相平 衡所需的分子動力學步。 表 3.2 不同分子數的計算時間 tab.3.2 the computer time of variable molecular number 分子數 n 10000 步所需時間(分鐘) 氣-液相平衡時間(分子動力學步) 1000 31 300000 2000 120 200000 3000 243 10000 4000 483 50000 5000 792 30000 從表 3.1、 表 3.2 可以看出： 計算同樣的分子動力學步(10000 步)， 分子數不同， 計算耗時也大相徑庭，分子數從 1000 到 5000，增加 5 倍，計算耗時則從 31 分鐘 到 792 分鐘，增加近 25 倍；對于達到氣-液相平衡所需的分子動力學步，則隨分子 數的增加，明顯減少，從 1000 個所需的 30 萬步減少到 5000 個分子數所需的 3 萬 步。綜合以上因素，本文確定：模擬分子數 n 的選取標準為：1000r * 時并不能忽略不計時而采取的一種校正方法。 長程校正的主要假設是徑向分布函數 rg(rr *)=1， 它對均勻相的校正非常適 用，對非均勻相，由于徑向分布函數 rg(rr *)=1 條件并不成立(我們將在第四章 討論液體、氣體、氣-液共存體的徑向分布函數時，專題說明)。因此，對氣-液界 面的分子動力學模擬采用長程校正方法效果并不明顯， m.mecke74等人對長程校正 方法應用于氣-液界面的分子動力學模擬的效果進行了研究，表 3.3 列出了 m.mecke 等人對 l-j 流體氬在不同溫度、 不同截斷半徑條件下， 采用長程校正(lrc) 與不采用長程校正(without lrc)的對比。 表 3.3 采用長程校正(lrc)與不采用長程校正(without lrc)的比較 tab.3.3 the comparison of simulation results of lrc and without lrc t r * lrc * l * v 0.70 2.5 no lrc 0.7861 0.00784 0.70 2.5 lrc 0.8154 0.00484 0.70 5.0 no lrc 0.8360 0.00194 0.70 5.0 lrc 0.8375 0.00201 0.85 2.5 no lrc 0.6956 0.03025 0.85 2.5 lrc 0.7391 0.01653 0.85 5.0 no lrc 0.7659 0.01041 0.85 5.0 lrc 0.7698 0.00977 分析表 3.3 可以發(fā)現(xiàn)： 當截斷半徑較小時， 采用長程校正， 均勻相中飽和氣體、 飽和液體密度以及非均勻相中表面張力有較大的改進，但與實驗值相比，仍然存 在較大的誤差；當截斷半徑較大時(通常r *5 后)，校正前后均勻相中飽和氣體、 3 氣-液界面微觀特性的分子動力學模擬研究 47 飽和液體密度以及非均勻相中表面張力改進并不明顯，反而由于增加了計算工作 量，采用長程校正方法的計算機用時大幅度增加。 3.1.3 氣-液界面的微觀特性 氣-液界面的微觀特性采用氣-液界面層厚度和氣-液分子數來定義，氣-液界面 層厚度和氣-液界面層內的分子通過下列兩種法則來確定： 10-90法則：氣-液界面層厚度根據如下公式確定： vlv %10 (3.10) vlv %90 (3.11) 即：將飽和氣體和飽和液體的過渡區(qū)域分成若干個等密度面，氣-液界面層由 介于等密度面 vlv %10和 vlv %90之間的區(qū)域構成，這 兩個等密度面之間的距離即是氣-液界面層的厚度。 由于一般情況下， 飽和氣體密度 * v 遠遠小于飽和液體密度 * l ， 所以通常也可 以認為氣-液界面層厚度是由 10%飽和氣體密度和 90%飽和液體密度之間的區(qū)域。 (25)/1.2rm法則80：(25)/1.2rm法則規(guī)定以 1.2 倍液相分子間的平均距 離為標準，對體系中的每個分子所處的狀態(tài)進行判斷，如果某個分子在該標準距 離內的其它分子數小于 2，則視為氣相分子，大于 5，則視為液相分子，在 25 之間則視為氣-液界面層上的分子。 通過本節(jié)的討論，綜合考慮各種因素。在本文的模擬中，均不考慮長程校正 對模擬結果的影響，除特別申明外，均將模擬盒沿 z 方向等分為 ns=100 個切片， 模擬分子數為 2048。由于體系達到平衡需要很長的模擬時間，為了得到較為真實 的平衡體系，本文取時間步長 t*=0.001，達平衡態(tài)的模擬長度為 10 萬分子動力學 步，然后以 12 萬分子動力學步進行統(tǒng)計。 3.2 新的氣-液界面微觀特性模擬算法變截斷半徑算法 目前氣-液界面的分子動力學模擬研究大多集中在一些宏觀參數，如：氣-液兩 相系統(tǒng)中飽和蒸汽壓力的計算、均勻相的飽和氣體和飽和液體密度的計算以及非 均勻相的表面張力的計算上，而對氣-液界面層的微觀結構如：界面層厚度、界面 層中的分子數，界面層中的分子漲落現(xiàn)象以及分子漲落現(xiàn)象對表面張力的影響研 究沒有給予足夠的重視，以至分子動力學模擬產生的大量的微觀信息被忽視，這 恰恰有悖于運用分子動力學模擬氣-液界面現(xiàn)象的初衷。 在本節(jié)里，本文圍繞氣-液界面的微觀特性展開研究：充分利用分子動力學模 擬過程中每個分子動力學步產生的微觀信息，并以此作為下一步計算的依據，以 此方法研究影響氣-液界面宏觀及微觀特性的主要因素，開發(fā)一種能真實反映氣- 重慶大學博士學位論文 48 液界面宏觀及微觀特性的分子動力學模擬算法。 3.2.1 傳統(tǒng)的模擬方法 在分子動力學模擬中，截斷半徑的選取對流體的模擬結果的影響非常重要， 從理論上講，截斷半徑越大，計算精度越高；但實際上：截斷半徑的選取受模擬 盒尺寸的限制，同時，截斷半徑越大，參與作用的分子數越多，計算速度越慢， 計算處理誤差越大。傳統(tǒng)的模擬方法在模擬過程中往往根據研究對象選取一個固 定截斷半徑，而不同的作者所推薦的截斷半徑不同且十分分散。由于氣-液相平衡 體系中，均勻相、非均勻相的分布特性不同，傳統(tǒng)的方法中統(tǒng)一的截斷半徑不可 能同時滿足均勻相中飽和氣體、飽和液體密度以及非均勻相中表面張力的計算精 度要求。下面，我們首先采用傳統(tǒng)的模擬方法-即固定截斷半徑對氣-液界面現(xiàn)象進 行模擬。 為了準確地反映截斷半徑對均勻相飽和氣體密度、飽和液體密度以及非均勻 相表面張力的影響， 本文首先采用傳統(tǒng)模擬方法， 以 l-j 流體 ar分子(參數見表 2.1) 為模擬對象，初始參數設置為：溫度 t=0.753，密度 * l =0.821，分子數為 2048， 取截斷半徑r *=2.5、3.0、3.5、4.0、4.5、5.5、6.5 等七種情況進行了分子動力學模 擬。表 3.4 表示了飽和氣體密度和飽和液體的密度以及表面張力模擬結果。 表 3.4 飽和氣體密度和飽和液體的密度以及表面張力模擬結果 tab.3.4 the simulation results of the density of saturated vapor/liquid and the surface tension r * * l * v * 2.5 0.7840 0.01144 0.5396 3.0 0.8125 0.00663 0.6183 3.5 0.8251 0.00460 0.7167 4.0 0.8345 0.00415 0.7683 4.5 0.8391 0.00371 0.8174 5.5 0.8416 0.00337 0.8393 6.5 0.8446 0.00330 0.8482 圖 3.2 為截斷半徑r *=2.5、3.5、4.5、5.5、6.5 條件下的氣-液態(tài)密度分布圖： 圖 3.3 是飽和液體密度和飽和氣體的密度及表面張力的模擬結果隨截斷半徑 3 氣-液界面微觀特性的分子動力學模擬研究 49 圖 3.3 汽相密度、液相密度、表面張力隨截斷半徑的變化 fig.3.3 density of saturated vapor/liquid and the surface tension vs truncated radii 的變化示意圖 圖3.2 不同截斷半徑下的氣-液平衡態(tài)密度分布 fig.3.2 density distributions of liquid-vapor equilibrium system vs truncated radii 工質：氬 分子數：n=2048 z* * 234567 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 md data r * l * v * * 重慶大學博士學位論文 50 從表 3.4 和圖 3.2、圖 3.3 可以看出： 在均勻體系中，飽和液體的密度 l 隨截斷半徑r *的增大而增大；飽和汽 體的密度 v 則隨r *的增大而減小。當截斷半徑較小時，兩者的模擬結果均與實驗 值相差較大，在r *=3.54.0時，模擬結果與實驗值(實驗值： l =0.821， v =0.00445)十分接近，飽和液體的密度誤差在 1%以內，飽和汽體的密度誤差在 4%以內。 在非均勻體系中，表面張力 的數值隨截斷半徑r *的增加而增加，并逐 漸向其實驗值逼近；當截斷半徑r *取值較小時，模擬結果的誤差較大；在截斷半 徑r *=3.54.0處，表面張力的計算誤差仍在 814%左右；在 r *=5.5處，表面 張力的模擬結果(=0.8393)與實驗值(實驗值：=0.8353)十分接近，但此處飽和汽 體的密度的計算誤差又明顯加大，達到約 25%。 不存在一個能同時滿足均勻相飽和氣體密度、 飽和液體密度以及非均勻相 表面張力特性計算的最佳的截斷半徑；對于均勻相、非均勻相各自的熱物理特性 計算，需取不同的截斷半徑，滿足非均勻相特性計算要求的截斷半徑應大于滿足 均勻相特性計算要求的截斷半徑。 分析圖 3.2 可以看出：勢能截斷半徑的選取還直接影響氣-液界面層的厚 度，根據按照10-90法則定義的氣-液界面區(qū)，則氣-液界面層的厚度隨勢能截斷 半徑的增加而減少，氣-液界面層的分子數也隨勢能截斷半徑的增加而減少。 3.2.2 變截斷半徑模擬方法 分析上面的模擬結果，我們認為：造成均勻相、非均勻相所要求的截斷半徑 不相同的主要原因是因為氣-液界面層內分子的不均勻分布造成。一般說來，運用 勢函數時， 我們能對分子i所受的來自于它截斷半徑內的所有分子對它的作用力作 準確計算。當截斷半徑較小時：在均勻相中，來自于截斷半徑以外的所有分子對 這個分子的附加作用力被平衡，故對它的運動軌跡沒有影響；而在非均勻相中， 自于截斷半徑以外的附加作用力不能被平衡，當求解運動方程時，這個附加作用 力必須考慮進去，所以，計算非均勻相的特性時，需要較大的截斷半徑。 a.lotfi 等早就注意到傳統(tǒng)的模擬方法的缺陷，第一次將長程修正(lrc： long-range corrections)的概念引入到氣-液界面特性的分子動力學模擬研究中； m.mecke 等74則在傳統(tǒng)的模擬結果的基礎上附加一個表面張力校正值 tail的方法 來計算表面張力； andrij, t.等3從勢函數截斷半徑和作用力函數截斷半徑不同步出 發(fā)， 提出 st(spherically truncated)勢函數模型和 sts(spherically truncated and shifted) 勢函數模型的概念，指出在進行傳統(tǒng)的模擬時，必須在截斷半徑處加上一個附加 作用力。 本文提出一種新的計算方法-變截斷半徑的計算方法，該方法克服了傳統(tǒng)模擬 3 氣-液界面微觀特性的分子動力學模擬研究 51 方法只重視結果統(tǒng)計，而忽視計算過程的缺點，充分利用模擬過程中已產生的分 子信息，定義10-90法則規(guī)定的區(qū)域為氣-液界面區(qū)，根據分子所處的位置改變截 斷半徑，當分子處于均勻相時，截斷半徑取一個固定的常數(通常取值較小)，當分 子處于非均勻相中，截斷半徑取分子到液相區(qū)邊界或汽相區(qū)邊界中距離較大者的 值，具體算法為： 1)根據模擬盒中分子分布的微觀信息，確定氣-液界面層的起始位置 1 z和終了 位置 v z； 2)根據分子位置z確定分子所處的狀態(tài)： 如果 v zz 或者 1 zz ，則分子處于均勻相-飽和氣相或飽和液相中； 如果 v zzz 1 ，則分子處于非均勻相-氣-液共存相中。 3)根據分子所處狀態(tài)，確定截斷半徑： 如果分子處于飽和氣相或飽和液相中，取定截斷半徑； 如果分子處于氣-液共存相中： (a) 計算zz 1 和 v zz ，比較zz 1 和 v zz 的大小； (b) 取截斷半徑為zz 1 和 v zz 中較大者。 (在上面的計算、比較中，由于氣-液界面層是對稱的，各項的數值均取正值) 在以上的計算中，由于氣-液界面區(qū)的邊界( 1 z和 v z)是變化的，同時分子所處 的位置 z 也是變化的，因而實現(xiàn)了模擬過程中的變截斷半徑算法。 對前述的模擬體系，本文采用變截斷半徑算法，在均勻相中截斷半徑分別取 r *=2.5、3.0、3.5、4.0 四種情況進行了模擬，飽和氣體密度和飽和液體