基于粗粒化分子動力學的自由水與水化硅酸鈣孔隙水凍結模擬

張 浩, 田 霞, 顧 鑫, 章 青

(河海大學 力學與材料學院,南京 211100)

1 引 言

以混凝土為代表的水泥基材料廣泛應用于各種工程中。在周期性變溫荷載作用下,水泥基材料孔隙內(nèi)的水將發(fā)生相變,導致水泥基材料損傷開裂、力學性能劣化并影響材料的服役性能。混凝土的凍融破壞在中國北方地區(qū)非常普遍[1],嚴重威脅混凝土結構的安全。因此,深入開展水泥基材料的凍融破壞機制研究[2],對于提升混凝土材料的抗凍性和耐久性,保障混凝土結構的安全具有重要的現(xiàn)實意義[3]。

水泥水化后生成水化硅酸鈣C-S-H(Calcium Silicate Hydrate)、氫氧化鈣和鈣礬石等產(chǎn)物,其中C-S-H占比60%～70%,極大地影響著水泥基材料的微觀結構特征和宏觀力學性能,故本文選用C-S-H表征水泥基材料的孔壁。Zhang等[4]建立了基于孔隙率變化的混凝土損傷模型,定量描述了宏觀力學性能的退化規(guī)律,提出凍融循環(huán)作用下混凝土動力增長系數(shù)的經(jīng)驗關系?；炷猎诮?jīng)歷凍融循環(huán)后,距離水面越近的部位孔隙率越大,其中小孔含量較少而大孔含量明顯增大,即距離水面越近凍融損傷越嚴重[5]。水泥基材料的凍融破壞主要是其孔隙中水的凍結膨脹產(chǎn)生的壓力所致[6],其凍融行為與受限水的相變特性直接相關,而受限水的相變特性遠比自由水復雜[7,8]。Foroutan等[9]使用分子動力學方法模擬了石墨烯納米孔隙中水的結冰過程,得到石墨烯納米孔隙通道中水的冰點為240 K,比自由水的冰點低29 K;Jiang等[10]通過分子動力學模擬發(fā)現(xiàn),6 ?納米硅孔隙中的液態(tài)水在2.5 GPa壓強和250 K溫度下,自發(fā)凍結成單層冰;8 ?納米硅孔隙中的液態(tài)水在4.0 GPa壓強和300 K溫度下,自發(fā)凍結成雙層冰。水泥基材料的凍融破壞研究表明,水泥孔隙分布對其損傷的形成十分重要。水泥基材料孔隙水凍結時的體積膨脹是導致水泥基材料凍融破壞的重要因素,研究表明受限水相對自由水冰點更低,更加難以結冰,但目前的研究中鮮有涉及水泥基材料孔隙水的凍結模擬。本文使用分子動力學方法圍繞不同孔徑的C-S-H孔隙水展開研究,模擬C-S-H孔隙水的凍結過程,探究不同孔隙對水泥基材料凍融破壞的危害程度。

雖然基于全原子分子動力學方法進行微觀尺度模擬精度較高,但對計算規(guī)模有較大限制,且很難模擬孔隙水的凍結過程,粗粒化分子動力學方法能夠有效地提升計算效率,可望進行孔隙水的凍結模擬。Jennings[11]提出了C-S-H球形膠體顆粒模型,采用Mie勢函數(shù)表示該粒子間的相互作用勢,并給出了從基本結構封裝形成更大尺度結構的方法。Ioannidou等[12]使用上述C-S-H球形膠體顆粒模型,計算得到了不同初始水灰比和密度下C-S-H樣品的局部壓力,分析了水灰比和密度對其力學性能的影響。Marrink等[13]提出脂質(zhì)的粗粒化模型,其中水的粗?；Ｗ邮菍⑺膫€水分子視為一個粗粒化水珠,即P4粒子,魏來[14]研究了該粗?；Ｗ觿莺瘮?shù)對結冰溫度的影響。從目前的研究來看,國內(nèi)外學者雖然分別建立了C-S-H和水的粗?；Ｐ?并圍繞相關問題開展研究,但鮮有同時建立這兩種粒子的粗?；Ｐ?以研究水泥基材料的凍融相關問題。

本文采用P4粒子和C-S-H膠體粒子,建立兩種粒子之間的耦合勢函數(shù),采用洛倫茲-貝特洛平均法則率定勢函數(shù)參數(shù)。基于上述兩種粗?；Ｗ?構建C-S-H孔隙及孔隙水的計算模型,模擬自由水與C-S-H孔隙水的結冰過程,研究C-S-H孔隙水的冰點與孔徑關系,并探究孔隙水結冰的微觀特性。

2 模型構建

水的粗?；Ｐ褪菍⑺膫€水分子合并為一個粗?；甈4粒子,則體系的粒子數(shù)量為全原子模型的1/4。粗?；疨4粒子間的相互作用勢采用LJ(Lennard-Jones)勢,即

(1)

式中r為粒子間距離,ε1為P4粒子間的勢阱深度,反映粒子間相互吸引作用的強弱,σ1為U1(r)=0時P4粒子間的距離。

將C-S-H粗?；汕蛐文z體顆粒,粒子間采用Mie勢函數(shù)[12],即

(2)

式中ε2為C-S-H膠體顆粒間的勢阱深度,σ2為U2(r)=0時C-S-H膠體顆粒間的距離。

將上述兩種勢函數(shù)進行加權平均,得到水的粗粒化P4粒子和C-S-H的粗?；Ｗ娱g的相互作用勢函數(shù),可表示為

Ucou(r)=λU1(r)+(1-λ)U2(r)

(3)

式中λ和1-λ分別為LJ勢函數(shù)和Mie勢函數(shù)的權函數(shù),λ=1,式(3)退化為P4粒子之間的LJ勢函數(shù);λ=0,式(3)退化為C-S-H膠體顆粒之間的Mie勢函數(shù)。

當λ=0.5時(后有兩種粒子相互作用準確性的驗證),將LJ和Mie勢函數(shù)表達式代入式(3),并將粒子間的勢阱深度和距離統(tǒng)一采用ε3和σ3表征,則有

(4)

式中ε3為P4粒子和粗粒化C-S-H粒子間的勢阱深度,σ3為U3(r)=0時兩種粒子間的距離。

兩種粒子之間的勢函數(shù)參數(shù)由洛倫茲-貝特洛法則確定,即

(5)

按上述方式形成水的粗?；疨4粒子和C-S-H的粗?；Ｗ娱g的相互作用勢函數(shù)后,便可構建C-S-H孔隙水的計算模型,離散后C-S-H孔隙通道的x-z截面如圖1所示。

圖1 C-S-H孔隙模型的x-z截面

3 自由水的結冰模擬

3.1 結冰判斷準則

本文采用鍵序參數(shù)描述水的結冰程度。鍵序參數(shù)[15]為每個粒子鍵方向順序的參數(shù),其取值范圍為0～1,0表示完全無序,1表示有序,其表達式為

(6)

(7)

冰是由水分子有序排列形成的結晶,其有序程度高于液態(tài)水。因此,可計算粗?；Ｗ拥逆I序參數(shù)均值,作為表征結冰的參數(shù)。

3.2 自由水凍結過程分析

首先,基于水的粗?；Ｗ幽Ｐ?計算粗?；幕疚锢硖匦?包括密度、表面張力與擴散系數(shù)),列入表1。由表1可知,數(shù)值模擬結果與理論值吻合較好。

表1 水的物理特性

為模擬自由水結冰的演化過程并分析自由水結冰過程中的變化特征,考慮邊長為10 nm的立方體模擬體系,隨機生成5000個P4粒子,采用等溫等壓系綜,溫度分別設置為270 K和275 K,壓強為標準大氣壓,為得到凍結過程中結冰演化過程細節(jié),模擬的時間步長選為較小的1 fs,截斷半徑為12 ?。參考魏來[14]判斷水是否結冰的方法,模擬結束后,統(tǒng)計等溫等壓系綜軌跡文件的最后一幀,觀察是否有晶體化現(xiàn)象產(chǎn)生,判斷水分子是否凍結。以此為依據(jù)選用0.44作為臨界值(未結冰的P4粒子鍵序參數(shù)在0.42以下,結冰后P4粒子鍵序參數(shù)在0.46以上):當鍵序參數(shù)小于0.44時,認為水處于液態(tài);當鍵序參數(shù)大于0.44時,認為水開始結冰。選用鍵序參數(shù)作為判斷結冰標準更加合理,且可以判斷局部位置是否結冰。圖2分別給出了270 K溫度下,7.7 ns,7.8 ns,7.9 ns和8.0 ns時x-z截面的鍵序參數(shù)云圖,為了更清楚地顯示自由水結冰的演化過程,將小于0.44的鍵序參數(shù)均設置為0。

從圖2可以看出,模擬時間小于7.7 ns時,計算區(qū)域的小部分粒子呈有序結構,大部分粒子呈無序結構;隨著模擬時間的增加,具有有序結構的區(qū)域迅速擴大,且呈規(guī)則排列;當模擬時間為8 ns時,計算區(qū)域基本均呈現(xiàn)有序結構,表明自由水已凍結成冰。在結冰過程中,水的勢能(1 Kc=4185.85 J)與流動性將發(fā)生變化,水的流動性可通過均方位移體現(xiàn),圖3和圖4分別給出了計算得到的勢能與均方位移曲線。

圖3 不同溫度下勢能隨時間的變化曲線

圖4 不同溫度下均方位移隨時間的變化曲線

從圖3和圖4可以看出,當系統(tǒng)溫度為275 K時,即不會發(fā)生凍結的情況下,水分子的勢能基本不發(fā)生變化(弛豫過程除外),均方位移與時間呈線性關系,符合水的基本特性。當系統(tǒng)溫度為270 K時,自由水將發(fā)生凍結,在凍結時,水分子的勢能將迅速降低,凍結后水分子的勢能和均方位移幾乎不再發(fā)生變化,表明體系中水已經(jīng)結冰,由液態(tài)變?yōu)楣虘B(tài),不再具有流動性。

模擬結果表明,粗粒化粒子凍結溫度與理論值273.15 K比較接近。此外,計算結果呈現(xiàn)了水凍結時的放熱過程以及從液態(tài)到固態(tài)的轉變,與Matsumoto凍結模擬的現(xiàn)象一致[16],表明粗粒化水P4用于凍結模擬的可行性。

4 C-S-H孔隙水的結冰模擬

4.1 C-S-H與水的相互作用驗證

自由水的冰點主要由水自身的性質(zhì)決定(即粗粒化P4粒子勢函數(shù)參數(shù)),而C-S-H孔隙水冰點受多種因素影響,除了水自身的性質(zhì)外,研究表明孔隙的孔徑[17]和壓強影響著孔隙水的冰點[18]。本文主要采用Lammps中fix NPT命令實現(xiàn)控溫控壓,保持壓強不變,研究C-S-H孔徑對該體系孔隙水冰點的影響。

在建立C-S-H孔隙結構計算模型時,上下區(qū)域由隨機生成的粗?；疌-S-H粒子組成,并設置粒子不發(fā)生移動,中間區(qū)域為隨機生成的粗粒化水,中間區(qū)域的高表示C-S-H的孔徑。設置計算體系為等溫等壓系綜,壓強為標準大氣壓,時間步長選用10 fs,模擬時三個方向均采用周期性邊界條件,粒子間的相互作用勢截斷半徑為12 ?。

考慮微孔半徑為5 nm的C-S-H孔隙結構,孔隙水在235 K結冰時,圖5和圖6分別給出該溫度下模擬得到的C-S-H孔隙結構x-z截面和y-z截面的粒子位置圖像。

圖5 235 K粒子位置圖像(x-z截面)

圖6 235 K粒子位置圖像(x-y截面)

圖7給出了相應的孔隙水結冰云圖。將圖7與圖2所示的自由水結冰云圖相比,可以看出受限水的結冰特征不同于自由水。對于受限水,在靠近C-S-H孔隙壁面的藍色區(qū)域,C-S-H吸附的水分子對應的鍵序參數(shù)小于0.44,呈現(xiàn)無序排列,未達到結冰狀態(tài);而在中間區(qū)域呈現(xiàn)有序排列,且鍵序參數(shù)大于0.44,達到結冰狀態(tài);而自由水的結冰區(qū)域覆蓋了x-z截面的絕大部分區(qū)域。

圖7 C-S-H孔隙x-z截面水結冰云圖

本文模擬結果表明,孔隙壁上會有一層吸附水膜,在極低的溫度下也不會結冰。這是因為靠近孔壁處,C-S-H與水的相互作用較大,這對水的運動有一定的束縛,從而導致C-S-H孔壁表面的水難以結冰,形成吸附水膜。C-S-H孔隙的區(qū)域2nm～12nm(紅線區(qū)域內(nèi)),模擬所得吸附水膜厚度約為0.5nm～0.8nm,與文獻[19]描述該溫度下吸附水膜厚度為0.63nm的情形相吻合。這表明該勢函數(shù)可以較為準確地模擬C-S-H與水的相互作用,而C-S-H對孔隙水冰點的影響主要考慮其表面張力,因此該模型適用于C-S-H后續(xù)計算。

4.2 C-S-H孔隙水密度分布

為探究C-S-H孔隙水的密度分布特征[20],將研究對象沿x和z方向均勻離散為邊長為δ的正方形,計算區(qū)域y方向長度作為高,形成有限個正四棱柱。統(tǒng)計每一幀正四棱柱內(nèi)粗粒化水分子的個數(shù),對時間平均,得出正四棱柱的數(shù)密度,再轉換為質(zhì)量密度。密度計算公式為

(8)

根據(jù)式(8)求出水的密度,可以得到C-S-H孔隙結構x-z截面水的密度云圖,如圖8所示。從圖5、圖6與圖8可以看出,在靠近C-S-H孔隙結構的邊界,受C-S-H表面張力的影響,水分子吸附在C-S-H粒子上,水分子密度較大,為1.5g/cm3～2g/cm3;在C-S-H孔隙結構的中間區(qū)域,水分子密度為1g/cm3左右,水分子呈現(xiàn)有序排列,已經(jīng)凍結。這表明水泥基材料孔隙中,C-S-H對水的影響隨著距離逐漸減小,C-S-H對孔壁表面處的水影響最大,對孔隙中間的水基本沒有影響。

4.3 不同孔徑中孔隙水的冰點計算

由于表面張力的作用,在水泥基材料不同孔徑的孔隙中,水的結冰溫度(冰點)有較大的差異。現(xiàn)有研究表明,孔徑越小,表面張力越大,冰點越低,可表示為[17]

圖8 C-S-H孔隙x-z截面水的密度云圖

(9)

Williams[21]實驗得出的-1 ℃下能結冰的最小半徑為100nm左右,與上述理論值的124.4nm較為接近;謝超等[22]推導了微孔半徑與冰點的關系,并參照該微孔半徑與冰點關系進行了實驗,該關系中-3 ℃下能結冰的最小半徑為41.4nm,與理論值41.78nm較為接近。兩者均與熱力學理論比較吻合,而更小尺度孔徑的冰點難以通過實驗獲得。

為研究微納米尺度下C-S-H孔徑對冰點的影響,分別設微孔半徑為5nm,10nm,15nm,20nm和25nm,受限水的模型上下邊界為3nm的C-S-H,中間區(qū)域生成的粗?；肿訑?shù)量依次為8000,16000,24000,32000和40000。計算不同溫度下的鍵序參數(shù),列入表2。以5nm微孔半徑C-S-H孔隙水為例,當溫度在242K及以下時,鍵序參數(shù)變化不大,且有序程度較好,表明孔隙水處于結冰狀態(tài);當溫度在242K以上時,鍵序參數(shù)迅速減小,有序程度較差,表明孔隙水處于自由狀態(tài)。據(jù)此,可以認為當鍵序參數(shù)在1K內(nèi)發(fā)生突變,則此時的溫度為孔隙水的冰點。

依據(jù)鍵序參數(shù)的結果可得到孔徑為5nm,10nm,15nm,20nm和25nm時,水的冰點依次為242K,251K,259K,262K和263K,數(shù)值模擬結果與式(9)熱力學理論結果對照如圖9所示。

表2 鍵序參數(shù)均值與溫度關系

表3 誤差分析

圖9 微孔半徑與冰點關系曲線

從圖9可以看出,模擬結果與理論值吻合較好,計算得到的不同孔徑的冰點略低于理論估算值,誤差分別為1.6%,3.5%,1.9%,1.8%和1.8%,最大誤差為3.5%。本文采用控溫控壓系綜保持模擬系統(tǒng)壓強不變,只考慮C-S-H孔徑對受限水冰點的影響,且選用的勢函數(shù)參數(shù)較為準確地表現(xiàn)了C-S-H與水的相互作用,故模擬結果與熱力學理論值相差較小。此外,研究發(fā)現(xiàn),當孔徑較小時,受限水的冰點波動較大;當孔徑較大時,冰點的變化緩慢。這表明水泥基材料中孔隙孔徑越小,其冰點越低,結冰難度會劇烈增加。因此,水泥基材料的孔隙大小對凍融破壞的影響巨大。由模擬所得的微孔半徑與冰點的關系可知,微孔半徑在5nm以下的孔在-31 ℃以上不會結冰,這意味著在通常的凍融環(huán)境條件下,孔徑在10nm以下的孔對水泥基材料基本沒有危害;微孔半徑在5nm～25nm的孔對應的孔隙水冰點在-31 ℃～-10 ℃,表明孔徑在10nm～50nm范圍內(nèi)的孔具有一定的危害;當微孔半徑在25nm以上時,孔隙水的結冰溫度在-10 ℃～0 ℃,表明孔徑在50nm以上的孔隙受凍融的影響較為明顯,對水泥基材料的危害也較大。本文的計算結果與熱力學理論相互印證,得出了實驗無法測得較小孔微孔半徑與冰點的關系,為進一步了解水泥基材料孔隙在凍融循環(huán)過程中的危害程度提供了基礎。

5 結 論

本文選取水和C-S-H的粗粒化粒子作為基本粒子,使用兩種粗?；Ｗ幽M自由水與C-S-H孔隙水的結冰過程,得出以下結論。

(1)C-S-H表面張力影響著孔隙水的密度分布和凍結分布特征。由于C-S-H壁面的表面張力作用,靠近C-S-H壁面處的水吸附在C-S-H粒子上,水的密度較大,且水的運動受到一定的束縛,形成無法結冰的吸附層;而在C-S-H孔隙結構中間區(qū)域,受C-S-H表面張力影響較小,該部位孔隙水的密度與自由水基本相同。當C-S-H孔隙較小時,結冰溫度非常低且吸附層占據(jù)主要區(qū)域,孔隙中水無法凍結。因此,在水泥基材料凍融破壞中,孔徑小的凝膠孔難以結冰,對水泥基材料的破壞基本沒有影響。

(2) 水泥基材料孔隙的孔徑對孔隙水冰點有著重大的影響,進而影響內(nèi)部孔隙水對水泥基材料凍融破壞的危害程度。水泥基材料孔隙的孔徑越小,孔隙水的冰點越低,且結冰難度會劇烈增加;孔徑在50nm以上時,C-S-H孔隙水的冰點在-10 ℃～0 ℃,孔隙水容易凍結,此類孔隙在凍融過程中危害相對較大;孔徑在10nm～50nm,孔隙水的冰點在-31 ℃～-10 ℃,此孔徑范圍的孔隙在凍融過程中具有一定的危害性;孔徑在10nm以下,孔隙水的冰點在-31 ℃以下,在凍融過程中基本不會結冰,危害很小。