含裂紋碳納米管軸壓性能的分子動力學(xué)模擬*

韓強(qiáng) 王彩紅 姚小虎

(華南理工大學(xué) 土木與交通學(xué)院，廣東 廣州510640)

碳納米管(CNTS)［1］具有優(yōu)異的力學(xué)、電學(xué)、熱學(xué)等性質(zhì)，在實(shí)驗(yàn)和理論研究中已廣泛用于復(fù)合材料，能夠極大改善材料的性能. 碳納米管優(yōu)異的力學(xué)性質(zhì)與其結(jié)構(gòu)特性密切相關(guān)，然而在生產(chǎn)過程中不可避免地會出現(xiàn)各種缺陷，特別是Stone-Wales(SW)及空位缺陷［2］，這些缺陷的存在會削弱碳納米管的拉壓等性能.碳納米管在軸壓荷載作用下發(fā)生屈曲行為［3］，辛浩等［4］采用Morse 勢函數(shù)，研究了含單雙原子及SW 缺陷碳納米管的軸壓屈曲性能，認(rèn)為各類缺陷均會顯著降低納米管的屈曲性能.Hirai 等［5］研究了沿軸向分布的空位缺陷對碳納米管壓縮性能的影響，發(fā)現(xiàn)在缺陷處存在應(yīng)力集中現(xiàn)象，軸壓屈曲荷載明顯減小. Huq 等［6］采用分子力學(xué)模擬，研究了SW 缺陷對單層碳納米管軸壓行為的影響.Zhang 等［7］運(yùn)用分子動力學(xué)模擬，研究了SW 及空位缺陷對單層碳納米管軸壓性能的影響，結(jié)果表明，缺陷大大削弱了碳納米管的臨界屈曲荷載.Ang 等［8］采用Compass 勢，研究了含1 －3個空位缺陷的單層碳納米管的軸壓屈曲行為. Kulathunga 等［9］采用Compass 勢函數(shù)，研究了含對稱與非對稱空位缺陷單層碳納米管的軸壓行為. Chowdhury等［10］采用Brenner 勢函數(shù)，研究了完善及含缺陷碳納米管的拉壓扭轉(zhuǎn)等力學(xué)性能，認(rèn)為空位缺陷對碳納米管的影響高于SW 缺陷.Vijayaraghavan 等［11］采用Airebo 勢函數(shù)，研究了含缺陷單層碳納米管軸壓行為.Kinoshita 等［12］采用Airebo 勢函數(shù)，分析了含SW 缺陷碳納米管的軸壓性能. Eftekhari 等［13］采用Tersoff 勢函數(shù)，分別研究了含空位及SW 缺陷的單層與雙層碳納米管的軸壓性能，研究表明，SW 缺陷對碳納米管屈曲性能的影響較空位缺陷顯著. 可見，已有的工作多是對含單雙空位及SW 缺陷的研究，關(guān)于含裂紋碳納米管軸壓性能的研究工作很少.

文中運(yùn)用分子動力學(xué)方法，對含裂紋扶手椅型碳納米管的軸壓性能進(jìn)行了模擬分析，重點(diǎn)考慮了裂紋尺寸及管長變化對碳納米管軸壓屈曲性能的影響，詳細(xì)分析了含不同裂紋尺寸碳納米管及不同長度納米管的軸壓變形特性.

1 分子動力學(xué)模擬

1.1 碳納米管的物理模型

文中以扶手椅型(7，7)及(10，10)碳納米管為研究對象，模擬過程中，碳納米管的長度變化范圍為6 ～30 nm. 在碳納米管中心附近區(qū)域，沿著與軸向成30°及90°方向去掉一定數(shù)目的碳原子，構(gòu)建出與軸向成不同傾角θ 的裂紋，設(shè)定初始裂紋尺寸為CL，在模擬過程中分別取為0.000、0.710、0.994及1.562 nm，其中裂紋尺寸為0.000 nm 表示完善的碳納米管.計(jì)算模型結(jié)構(gòu)如圖1 所示.

圖1 長度為6 nm 的(7，7)碳納米管的結(jié)構(gòu)Fig.1 Structures of (7，7)carbon nanotube with a length of 6nm

1.2 分子動力學(xué)模擬方法

分子動力學(xué)計(jì)算的一個關(guān)鍵問題是原子勢函數(shù)的選取，文中選用Airebo 勢函數(shù)［14］，這是一種改進(jìn)的Rebo 勢函數(shù)，能夠較為真實(shí)地反映碳元素所構(gòu)成的固態(tài)材料的物理性質(zhì).

模擬過程中，控制X、Y、Z 方向?yàn)樽杂蛇吔鐥l件，沿著碳納米管軸向?yàn)榧虞d方向，固定一端碳原子，對另一端碳原子施加恒速壓縮位移載荷. 采用Nose-Hoover 熱浴法［15-16］進(jìn)行溫度調(diào)節(jié)，使系統(tǒng)保持在特定的溫度下，采用Velocity-Verlet 算法［17］進(jìn)行分子動力學(xué)計(jì)算.加載之前先對碳納米管進(jìn)行充分無約束弛豫，使整個系統(tǒng)處于能量最低的平衡狀態(tài).對弛豫過的碳納米管施加軸壓位移載荷，每步應(yīng)變量為0.000 5，然后進(jìn)行充分的弛豫，重復(fù)此加載、弛豫過程，直至碳納米管壓縮破壞.

2 模擬結(jié)果分析

2.1 裂紋尺寸對碳納米管軸壓性能的影響

為研究裂紋尺寸變化對碳納米管軸壓屈曲性能及變形特性的影響，選取長度為6 nm、初始裂紋尺寸分別為0.0、0.71、0.994 及1.562 nm 的(7，7)與(10，10)碳納米管進(jìn)行模擬分析，模擬過程中控制溫度為0.01 K，以避免熱激活的復(fù)雜影響，時間步長取1 fs.

2.1.1 垂直于軸線方向的裂紋

采用分子動力學(xué)方法，模擬含有垂直于軸向裂紋的(7，7)、(10，10)扶手椅型碳納米管在軸壓載荷作用下的屈曲及后屈曲變形行為.為了進(jìn)行對比分析，首先模擬了完善(7，7)管的軸壓屈曲變形性能，其納米管直徑約0.95 nm.

文中與辛浩［4］及Yakobson 等［3］所選取的勢函數(shù)、分子動力學(xué)計(jì)算方法及模型尺寸不完全相同，將他們的結(jié)果與文中的模擬結(jié)果進(jìn)行對比分析可知(如表1 所示)，文中得到的(7，7)管軸壓臨界屈曲應(yīng)變與他們的結(jié)果很接近，可見Airebo 勢函數(shù)可以較好地描述碳納米管中的C—C 鍵作用.

表1 完善(7，7)碳納米管在軸壓變形過程中關(guān)鍵點(diǎn)的應(yīng)變值Table 1 Critical strains of perfect (7，7)CNTs during axial compression

由圖2 可知，對于完善碳納米管，隨著壓縮進(jìn)行，其軸向應(yīng)力呈非線性變化.圖2 中a、b、c、d 四點(diǎn)與圖3(a)、3(b)、3(c)、3(d)一一對應(yīng)，當(dāng)軸壓應(yīng)變?yōu)?.1347(圖中d 點(diǎn))時，碳納米管徹底壓縮破壞.

圖2 完善(7，7)碳納米管的軸壓應(yīng)力－應(yīng)變曲線Fig.2 Stress-strain curve of perfect (7，7)carbon nanotube under axial compression

圖3 完善(7，7)碳納米管的軸壓變形Fig.3 Configurations of perfect (7，7)carbon nanotube under axial compression

圖4 θ=90°、裂紋尺寸為0.71nm 的(7，7)管的軸壓應(yīng)力－應(yīng)變曲線Fig.4 Stress-strain curve of (7，7)tube with the crack length of 0.71 nm under axial compression when θ=90°

圖5 θ=90°、裂紋尺寸為0.71 nm 的(7，7)管的軸壓變形Fig.5 Configurations of (7，7)tube with the crack length of 0.71 nm under axial compression when θ=90°

圖4 和5 給出了裂紋尺寸為0.71 nm 的(7，7)碳納米管在軸壓過程中的應(yīng)力－應(yīng)變曲線及變形形態(tài). 由圖5可知，在軸壓過程中，含裂紋(7，7)管的變形形態(tài)與完善(7，7)管有明顯區(qū)別，對于完善碳納米管(如圖3 所示)，隨著壓縮進(jìn)行，管壁某處發(fā)生凹陷后，凹陷會沿著軸向進(jìn)行傳播，管壁發(fā)生凹陷的地方可能會回彈.然而對于裂紋尺寸為0.71 nm的納米管，在整個軸壓過程中，裂紋處的凹陷沒有發(fā)生回彈，這是由于裂紋處缺少相應(yīng)的支撐作用.

裂紋附近區(qū)域受力不均勻，在壓縮至a 點(diǎn)時，裂紋區(qū)域出現(xiàn)局部凹陷變形.與完善(7，7)管相比，裂紋尺寸為0. 71 nm 時，其臨界屈曲強(qiáng)度下降了40.4%，當(dāng)壓縮至b 點(diǎn)時，與裂紋部位相對應(yīng)的另一側(cè)管壁出現(xiàn)凹陷. 對比完善管及含裂紋(7，7)管的應(yīng)力－應(yīng)變曲線可知，在b、c、d 處兩種管內(nèi)的應(yīng)力差別較小，可見裂紋尺寸為0.71 nm 時，納米管屈曲后軸向承載力的變化不大.碳納米管最終在裂紋處彎折破壞.

圖6、7 給出了裂紋尺寸為0.994 nm 的(7，7)管的軸壓性能. 由圖7 可知，裂紋尺寸為0. 994 與0.71 nm的(7，7)管的軸壓過程相似，在裂紋附近處均出現(xiàn)凹陷并最終在裂紋處壓縮破壞，在整個軸壓過程中，裂紋處的凹陷沒有發(fā)生回彈.由圖6可知，碳納米管在軸壓過程中荷載變化較為頻繁，在未達(dá)到臨界屈曲強(qiáng)度之前出現(xiàn)a、b 兩個波動點(diǎn)，兩點(diǎn)對應(yīng)的變形結(jié)構(gòu)表明含裂紋一側(cè)的管壁出現(xiàn)凹陷.壓縮至c 點(diǎn)時，管壁另一側(cè)也出現(xiàn)凹陷，此時納米管達(dá)到軸壓屈曲強(qiáng)度. 與完善管相比，其屈曲強(qiáng)度降低了45.4%.

圖6 θ =90°、裂紋尺寸為0. 994 nm 的(7，7)管的軸壓應(yīng)力－應(yīng)變曲線Fig.6 Stress-strain curve of (7，7)tube with the crack length of 0.994 nm under axial compression when θ=90°

圖7 θ=90°、裂紋尺寸為0.994 nm 的(7，7)管的軸壓變形Fig.7 Configurations of (7，7)tube with the crack length of 0.994 nm under axial compression when θ=90°

對比裂紋尺寸為0.994 與0.71 nm 的(7，7)管的變形可知，隨著裂紋尺寸增加，裂紋處的凹陷略有差別，對于裂紋尺寸較大的納米管，其凹陷處裂紋上下邊緣呈現(xiàn)明顯的錯位，這是由于裂紋處缺少了沿軸向的支撐作用，裂紋尺寸增大時，這種缺失越明顯.

圖8 給出了(7，7)、(10，10)管的臨界屈曲強(qiáng)度隨裂紋尺寸的變化，隨著裂紋尺寸增加，碳納米管的屈曲強(qiáng)度均減小. 在裂紋尺寸由0 增加到1.562 nm的過程中，(7，7)與(10，10)管的軸向臨界屈曲強(qiáng)度分別降低了57.05%及49.98%.

圖8 θ=90°時臨界屈曲強(qiáng)度隨裂紋尺寸的變化Fig.8 Critical buckling strength vs.crack length when θ=90°

2.1.2 與軸線方向成30°角的裂紋

由上節(jié)可知，碳納米管中部區(qū)域垂直于軸線的裂紋大大減弱了其軸向承載力，為了進(jìn)行更加全面的分析，本節(jié)模擬分析了裂紋傾角為30°、管長為6 nm，含不同裂紋尺寸的(7，7)、(10，10)扶手椅型碳納米管的軸壓屈曲性能.

圖9、10 給出了裂紋傾角為30°、尺寸為0.71 nm的(7，7)碳納米管的軸壓性能. 裂紋附近區(qū)域的管壁受力不均勻，碳納米管壓縮至a 點(diǎn)時，在裂紋附近的管壁首先發(fā)生凹陷變形，與完善管相比，其臨界屈曲強(qiáng)度降低了40.6%. 隨著荷載增加，管壁另一側(cè)也出現(xiàn)凹陷(對應(yīng)b 點(diǎn)形態(tài)圖)，其后的變形及對應(yīng)的軸向應(yīng)力及勢能變化與裂紋傾角為90°、尺寸為0.71nm 的(7，7)納米管相似，軸壓過程中相應(yīng)的軸向應(yīng)力小于裂紋傾角為90°的納米管.

裂紋傾角為90°的(7，7)納米管的凹陷出現(xiàn)在裂紋處，整個裂紋區(qū)域內(nèi)凹，而裂紋傾角為30°、尺寸為0.71 nm 的(7，7)納米管最開始出現(xiàn)的凹陷是在裂紋附近區(qū)域，而不是裂紋區(qū)域內(nèi)凹，在裂紋邊緣處出現(xiàn)扭曲現(xiàn)象. 對比二者的結(jié)構(gòu)可知，當(dāng)裂紋與軸向成30°角時，在軸向壓縮作用下，裂紋的上、下兩邊緣處有斜向剪切作用，這個剪切作用引發(fā)裂紋區(qū)域的扭曲，可見裂紋結(jié)構(gòu)差異會導(dǎo)致碳納米管變形形態(tài)及性能的不同.

圖9 θ=30°時裂紋尺寸為0. 71 nm 的(7，7)管的軸壓應(yīng)力－應(yīng)變曲線Fig.9 Stress-strain curve of (7，7)tube with the crack length of 0.71 nm under axial compression when θ=30°

圖10 θ=30°、裂紋尺寸為0.71nm 的(7，7)管的軸壓變形Fig.10 Configurations of (7，7)tube with the crack length of 0.71 nm under axial compression when θ=30°

裂紋傾角為30°、尺寸為0.994 nm 的(7，7)碳納米管的軸壓性能如圖11、12 所示.壓縮至a 點(diǎn)時，納米管達(dá)到其臨界屈曲強(qiáng)度，與完善管相比，其值降低了51%. 由圖可知，起初其變形形態(tài)及應(yīng)力變化趨勢與裂紋尺寸為0.71nm 的(7，7)管相似，在裂紋附近區(qū)域出現(xiàn)內(nèi)凹，裂紋上下邊緣出現(xiàn)扭曲現(xiàn)象.然而隨著載荷增加，可以明顯看到裂紋處出現(xiàn)劇烈的錯位及撕扯，其縱軸線發(fā)生明顯的彎折.

圖11 θ=30°、裂紋尺寸為0.994nm 的(7，7)管的軸壓應(yīng)力－應(yīng)變曲線Fig.11 Stress-strain curve of (7，7)tube with the crack length of 0.994 nm under axial compression when θ=30°

圖12 θ=30°、裂紋尺寸為0.994 nm 的(7，7)管的軸壓變形Fig.12 Configurations of (7，7)tube with the crack length of 0.994 nm under axial compression when θ=30°

圖13 θ=30°時臨界屈曲強(qiáng)度隨裂紋尺寸的變化Fig.13 Critical buckling strength vs the crack length when θ=30°

圖13 給出了裂紋傾角為30°的(7，7)、(10，10)碳納米管的臨界屈曲強(qiáng)度隨裂紋尺寸的變化，由圖可知，隨著裂紋尺寸增加，其臨界屈曲強(qiáng)度降低，在裂紋尺寸由0 增加到1.562 nm 的過程中，(7，7)與(10，10)管的臨界屈曲強(qiáng)度分別降低了58.45%及52.48%.

2.2 不同長度含裂紋碳納米管的軸壓屈曲

2.2.1 垂直于軸線方向的裂紋

對管長范圍為6 ～30 nm 的含垂直于軸線方向裂紋的(7，7)、(10，10)碳納米管進(jìn)行分子動力學(xué)模擬，分析其軸壓屈曲和后屈曲性能.

圖14、15 給出了長度分別為10 及20 nm、裂紋尺寸為0.994 nm 的(7，7)管在軸壓過程中的應(yīng)力－應(yīng)變曲線及變形形態(tài).與前面管長為6 nm 的管進(jìn)行對比分析可知，納米管的彎折同樣是首先出現(xiàn)在裂紋處.管長為10 及20 nm 的(7，7)管的軸壓屈曲性能的區(qū)別不大. 根據(jù)勢能變化及相應(yīng)的結(jié)構(gòu)形態(tài)(圖14(b)及圖15(b))可知，碳納米管的壓縮與壓桿的屈曲相似，在整體發(fā)生屈曲后并沒有出現(xiàn)局部的屈曲半波，隨著荷載增加，納米管多處出現(xiàn)彎折.

圖14 θ =90°、裂紋尺寸為0. 994 nm、管長為10 nm 的(7，7)管的軸壓Fig.14 Curves of stress-strain and potential energy of (7，7)tube with the crack length of 0.994 nm and the tube length of 10 nm under axial compression when θ=90°

圖15 θ=90°、裂紋尺寸為0.994nm、管長為20nm 的(7，7)管的軸壓應(yīng)力－應(yīng)變和勢能變化曲線Fig.15 Curves of stress-strain and potential energy of (7，7)tube with the crack length of 0.994 nm and the tube length of 20 nm under axial compression when θ=90°

由圖16 可知，對于含垂直軸向裂紋的(7，7)、(10，10)管，隨著管長增加，納米管的臨界屈曲荷載呈減小的趨勢，管長大于10 nm 后，兩種管徑的臨界屈曲荷載變化都不大.由圖16(a)可知，對于裂紋尺寸分別為0.994 和1.562 nm 的(7，7)納米管，隨著管長從6 增加到30nm，其臨界屈曲強(qiáng)度的減小量分別為32.8%和25.0%.由圖16(b)知，對于裂紋尺寸分別為0.994 和1.562 nm 的(10，10)納米管，隨著管長增加，其臨界強(qiáng)度分別減小了11.9%和13.4%.

2.2.2 與軸線方向成30°角的裂紋

本節(jié)模擬了裂紋傾角為30°、管長變化范圍為6 ～30 nm 的(7，7)碳納米管的軸壓變形行為.

圖17、18 給出了管長為10 及20 nm 的(7，7)管的軸壓性能.與管長為6 nm 的管進(jìn)行對比可知，納米管的凹陷同樣是首先出現(xiàn)在裂紋附近區(qū)域，裂紋上、下邊緣出現(xiàn)扭曲. 管長為10 及20 nm 的含裂紋(7，7)管的軸壓屈曲行為相似，由勢能變化及結(jié)構(gòu)形態(tài)可知，碳納米管的壓縮與壓桿的屈曲相似，在整體發(fā)生屈曲后并沒有出現(xiàn)局部的屈曲半波，隨著荷載增加，多處出現(xiàn)彎折，這一整體彎曲現(xiàn)象與裂紋傾角為90°、管長為10 及20 nm 的管相似，但相應(yīng)的裂紋附近區(qū)域的彎曲細(xì)節(jié)不同. 由圖17(a)與18(a)可知，對于裂紋傾角為30°的(7，7)管，隨著管長增加，其臨界屈曲強(qiáng)度有一定的降低，對于較長的碳納米管，其呈現(xiàn)壓桿屈曲的特性越明顯.

圖16 含裂紋管的臨界屈曲強(qiáng)度隨管長的變化Fig.16 Critical buckling strength of tube with crack vs. tube length

圖17 θ=30°、裂紋尺寸為0.994nm、管長為10nm 的(7，7)管的應(yīng)力-應(yīng)變和勢能變化曲線Fig.17 Curves of stress-strain and potential energy of (7，7)tube with the crack length of 0.994 nm and the tube length 10 nm under axial compression when θ=30°

圖18 θ=30°、裂紋尺寸為0.994nm、管長為20nm 的(7，7)管的應(yīng)力－應(yīng)變和勢能變化曲線Fig.18 Curves of stress-strain and potential energy of (7，7)tube with the crack length of 0.994 nm and the tube length of 20 nm under axial compression when θ=30°

3 結(jié)論

運(yùn)用分子動力學(xué)方法，研究了含不同裂紋的(7，7)、(10，10)碳納米管的軸壓屈曲行為，分別分析了裂紋尺寸及管長變化對扶手椅型納米管軸壓屈曲性能的影響，得到了相應(yīng)的應(yīng)力－應(yīng)變曲線及變形形態(tài)圖.結(jié)果表明，隨著裂紋尺寸增加，含裂紋碳納米管的臨界屈曲強(qiáng)度明顯降低.對于裂紋傾角為90°的納米管，在軸壓荷載作用下，裂紋區(qū)域內(nèi)凹，并呈現(xiàn)一定程度的錯位現(xiàn)象;對于裂紋傾角為30°的納米管，隨著荷載增加，裂紋附近區(qū)域出現(xiàn)凹陷，在裂紋的上下邊緣出現(xiàn)扭曲. 對于含兩種傾角裂紋的納米管，在整個軸壓過程中，裂紋附近區(qū)域的凹陷不再會發(fā)生回彈. 隨著管長增加，含裂紋碳納米管的軸向承載力同樣減小. 對于裂紋傾角為90°、裂紋尺寸為0.994 nm 的(7，7)和(10，10)碳納米管，在管長從6 增加到30 nm 的過程中，其臨界屈曲強(qiáng)度分別減小了32.8%及11.9%.另外由納米管相應(yīng)的變形形態(tài)可知，裂紋傾角變化也會影響納米管結(jié)構(gòu)的壓縮行為.

［1］Iijima S. Helical microtubules of graphitic carbon ［J］.Nature，1991，354(6348):56-58.

［2］Ebbesen T，Takad T. Topological and SP 3 defect structures in nanotubes［J］.Carbon，1995，33(7):973-978.

［3］Yakobson B I，Brabec C，Bernholc J. Nanomechanics of carbon tubes:instabilities beyond linear response ［J］.Physical Review Letters，1996，76(14):2511-2514.

［4］辛浩.石墨烯/碳納米管力學(xué)性能的研究［D］. 廣州:華南理工大學(xué)土木與交通學(xué)院，2010.

［5］Hirai Y，Nishimaki S，Mori H，et al. Molecular dynamics studies on mechanical properties of carbon nano tubes with pinhole defects ［J］. Japanese Journal of Applied Physics，2003，42(6S):4120.

［6］Huq A，Goh K，Zhou Z，et al. On defect interactions in axially loaded single-walled carbon nanotubes［J］. Journal of Applied Physics，2008，103(5):054306.

［7］Zhang Y，Xiang Y，Wang C.Buckling of defective carbon nanotubes ［J］. Journal of Applied Physics，2009，106(11):113503.

［8］Ang K，Kulathunga D，Reddy J. In buckling of defective carbon nanotube［C］∥Proceedings of the International Conference on Computing in Civil and Building Engineering.［S.l.］:Nottingham University Press，2010.

［9］Kulathunga D，Ang K，Reddy J.Molecular dynamics analysis on buckling of defective carbon nanotubes［J］.Journal of Physics:Condensed Matter，2010，22 (34):345301.

［10］Chowdhury S C，Haque B Z，Gillespie Jr J W，et al.Molecular simulations of pristine and defective carbon nanotubes under monotonic and combined loading［J］.Computational Materials Science，2012，65:133-143.

［11］Vijayaraghavan V，Wong C. Temperature，defect and size effect on the elastic properties of imperfectly straight carbon nanotubes by using molecular dynamics simulation［J］.Computational Materials Science，2013，71:184-191.

［12］Kinoshita Y，Kawachi M，Matsuura T，et al.Axial buckling behavior of wavy carbon nanotubes:A molecular mechanics study［J］.Physica E:Low-dimensional Systems and Nanostructures，2013，54:308-312.

［13］Eftekhari M，Mohammadi S，Khoei A R.Effect of defects on the local shell buckling and post-buckling behavior of single and multi-walled carbon nanotubes［J］.Computational Materials Science，2013，79:736-744.

［14］Stuart S J，Tutein A B，Harrison J A.A reactive potential for hydrocarbons with intermolecular interactions ［J］.The Journal of Chemical Physics，2000，112(14):6472-6486.

［15］Hoover W G. Canonical dynamics:equilibrium phasespace distributions ［J］. Physical Review A，1985，31(3):1695.

［16］Nosé S. A unified formulation of the constant temperature molecular dynamics methods ［J］. The Journal of Chemical Physics，1984，81(1):511-519.

［17］Swope W C，Andersen H C，Berens P H，et al. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules:application to small water clusters［J］. The Journal of Chemical Physics，1982，76(1):637-649.