濕度對水化硅酸鈣力學性能影響的分子動力學模擬

劉士達，李宗利,2，童濤濤，李云波，肖帥鵬

(1.西北農林科技大學水利與建筑工程學院，楊凌 712100；2.西北農林科技大學旱區農業水土工程教育部重點實驗室，楊凌 712100)

0 引 言

混凝土具有多尺度結構特征，不同尺度性能對混凝土宏觀性能影響不同?；炷敛牧系牧W性能主要決定于水泥水化產物，而水泥水化產物中最重要的成分是水化硅酸鈣(calcium silicate hydrate， C-S-H)凝膠體，占硬化水泥漿體體積的60%～70%，對水泥基材料的行為起著決定性的作用[1]。水化硅酸鈣是水泥基材料組成的基本單元，其特性是解釋混凝土宏觀性能的基礎。

水泥基材料經常服役于不同濕度的環境中，在水力梯度或自由吸水作用下，外界水進入使其處于不同的濕度狀態。水泥基材料中的水分以多種形式存在[1-2]，按照水的存在形式，可以分為自由水、毛細水和層間水。其中，層間水是C-S-H層間吸附的水分子，與自由水和毛細水一起隨濕度而變化。層間水對C-S-H力學性能的影響，一定程度上反映了濕度對水泥基材料的力學性能影響規律[3]。一方面，凝結硬化過程中層間水因水泥漿的初始含水量或水灰比大小在變化，從而影響水泥水合物的微觀結構。另一方面，使用過程中層間水對周圍環境和濕度也很敏感[4]，譬如，當浸入水中時，C-S-H凝膠吸水體積膨脹；當外界環境濕度降低時，C-S-H會脫去層間水，引起干燥收縮、斷裂等微觀結構的變化，從而影響C-S-H力學性能，而且不同受力狀態影響機理也不同，因此，有必要研究不同受力狀態下濕度對C-S-H力學性能的影響。

目前關于濕度對水泥基材料力學性能的研究大都基于宏觀或細觀尺度進行，微觀或分子尺度下濕度對水泥基材料力學性能影響已有少量研究成果。張偉等[5]采用反應力場的分子動力學方法研究了缺陷內水分子對C-S-H性能的影響。得到的結果是初始缺陷內的水不僅自身會形成氫鍵，也會與C-S-H基體形成氫鍵，略微提高C-S-H的抗拉強度和彈性模量。Hou等[6]采用分子動力學方法研究了不同濕度C-S-H在單軸拉伸和單軸壓縮下的力學性能。結果表明，濕度會降低C-S-H拉伸、壓縮時的力學性能，且C-S-H的抗壓強度遠大于抗拉強度。Zhou等[7]對飽和度為0%～100%的C-S-H進行了不同應變率的單軸拉伸分子動力學模擬。計算結果表明，C-S-H的動態拉伸力學性能隨著濕度的增加而降低。以上科研人員對不同濕度C-S-H的力學性能做了研究，然而僅僅針對的是C-S-H在受拉和受壓狀態下的力學性能。不同受力狀態下C-S-H表現出的力學特性不同，濕度對不同受力狀態下力學性能的影響程度也不同，需要進一步系統研究。

本文利用分子動力學，首先分析了不同濕度的C-S-H在結構和彈性參數上的差異，然后計算了不同濕度C-S-H在拉伸、壓縮、剪切不同受力狀態下的應力應變關系，對比分析了不同受力狀態下濕度的影響程度，拓寬濕度對混凝土力學性能影響研究范圍，系統揭示了納觀尺度下濕度對C-S-H力學性能的影響規律。

1 模擬方法

1.1 力場的選擇和模型的建立

力場是表示分子系統勢能變化的函數，選擇合理與否影響著仿真精度的高低[8]。水化硅酸鈣力場(CSHFF)是基于黏土力場(ClayFF)[9]開發的，主要通過改進ClayFF中Ca、Si和O的參數，使C-S-H的結構數據和高階性質更加精確[10]。該力場已被證明可以正確預測實際水泥水合物的基本結構特征和基本物理特征，并且在C-S-H力學性能的分子模擬中已被廣泛認可和應用[6,11-12]。因此，本研究選用CSHFF力場作為分子模擬的力場。通過文獻[10]可獲得CSHFF中描述Ca、Si、O、H原子間相互作用的參數。

本文基于Pellenq等[13]所提供的方法建立C-S-H模型。該方法基于Tobermorite模型刪除部分SiO2所建立的“真實”模型能夠很好地匹配實際C-S-H的各種參數，在C-S-H分子模擬中被廣泛應用[6,8,14]。已有研究證明該模型能夠很好地預測C-S-H的物理力學性能。模型建立的具體過程為：

(1)以Tobermorite 11?晶體作為C-S-H模型的初始結構，根據核磁共振(NMR)測試提供的Qn分布進行斷鏈[15]，即Q0=0%，Q1=66.01%，Q2=33.99%，平均鏈長(MCL=2(Q2/Q1+1)=3.03)與NMR測試結果保持一致，并參照Allen等[16]通過實驗所測得的C-S-H的平均分子式，鈣硅比采用1.7。然后在300 K、等溫等壓(NPT)系綜下弛豫，得到干燥的C-S-H。

(2)采用蒙特卡洛方法對干燥的C-S-H進行吸水，吸附過程在300 K、NPT系綜下完成，吸附步數設置為一億步。將壓強設置為飽和蒸汽壓的一半，當吸附達到平衡狀態時，得到濕度為50%的C-S-H；當吸附在飽和蒸汽壓下達到平衡狀態時，得到飽和狀態下的C-S-H。最終建立的模型在分子式、密度、Qn物種分布和平均鏈長等方面與實驗結果接近[15-16]。

(3)為了得到穩定統計結果，獲得可靠的實效模式[17]，將上述建立的不同濕度的C-S-H模型擴展成1×3×4的超胞。濕度為0%的C-S-H超胞的尺寸為2.0 nm×7.8 nm×8.0 nm，含有8 616個原子；濕度為50%的C-S-H超胞尺寸為2.1 nm×8.2 nm×8.4 nm，含有11 892個原子；濕度為100%的C-S-H超胞的尺寸為2.2 nm×8.6 nm×8.8 nm，含有15 168個原子。

圖1展示了不同濕度C-S-H的建立過程。

圖1 C-S-H模型的建立過程Fig.1 Establishment process of C-S-H model

1.2 分子動力學模擬

受目前常用分子模擬軟件分析功能限制，C-S-H彈性力學參數和不同受力狀態的模擬分別采用不同軟件進行。

C-S-H的彈性常數利用Materials Studio[18]軟件進行計算。計算過程中通過對C-S-H結構施加應變，分析結構變形得到C-S-H的彈性常數。C-S-H的體積模量、剪切模量根據彈性常數和柔性常數，采用Voigt-Reuss-Hill平均值法計算，楊氏模量根據柔性常數進行計算。

LAMMPS(大規模原子/分子并行模擬器)是一個經過良好測試和廣泛使用的開源經典分子動力學軟件[19]，在模擬單軸拉伸和單軸壓縮方面較為方便。對1.1節所建立的不同濕度C-S-H超胞模型進行充分弛豫，作為單軸拉伸、壓縮模擬的模型。拉壓沿平行和垂直于C-S-H的鏈狀方向進行，即y和z方向。模型三個方向均采用周期性邊界條件，模擬過程中對模型施加應變，應變速率保持0.08 ps-1不變。當沿一個方向進行拉伸或壓縮時，其他兩個方向的壓力保持為零，變形自由，體現泊松效應。垂直于加載方向的無壓力設置消除了對變形的人為約束，允許張力自由發展，不受任何約束[6]。整個系統的弛豫和模擬均在NPT系綜下進行，模擬時間步長設置為1.0 fs。

C-S-H的剪切模擬也采用Materials Studio軟件[18]完成。模擬剪切的模型采用1.1節所建立的C-S-H超胞。模型的三個方向均采用周期性邊界條件，模擬過程中對超胞在y-z平面上施加剪應變，應變速率為0.005 ps-1。在每次施加完成后，對超胞進行弛豫，弛豫時間為10 ps，并將充分弛豫后的模型作為下一次施加的初始模型。剪切和弛豫過程均采用NVT系綜，模擬時間步長設置為1.0 fs。每施加一次應變，都會得到一個相對應的剪應力，從而得到剪切時的應力和應變曲線。通過梯形積分規則計算應力應變曲線下的面積，得到拉伸、壓縮和剪切時的單位體積應變能密度，用應變能密度作為評價系統變形性能的間接指標。

2 結果與討論

2.1 不同濕度C-S-H結構及彈性常數的差異

圖2 不同濕度C-S-H超胞Si—O、Cas—O、 Caw—O徑向分布函數Fig.2 Radial distribution functions of Si—O, Cas—O and Caw—O in C-S-H supercell with different humidity

徑向分布函數(RDF)可以解釋為系統的區域密度與平均密度的比，反映微觀結構變化。對1.1節所建立的不同濕度的C-S-H超胞模型進行充分弛豫后，計算層間鈣與橋接氧Caw—O、層內鈣與橋接氧Cas—O、硅原子與橋接氧Si—O的RDF，如圖2所示。RDF的峰值所對應的位置表示兩種原子之間的鍵長。本文所建立的不同濕度C-S-H模型中硅原子與氧原子的鍵長為(0.161±0.01) nm，層內鈣與氧原子的鍵長為(0.231±0.02) nm，層間鈣與氧原子的鍵長為(0.237±0.02) nm。Richardson[20]的實驗數據顯示，Si—O平均鍵長為0.164 nm，Ca—O鍵長在0.22～0.29 nm；X光散射實驗[21]表明，Si—O平均鍵長為0.167 nm，Ca—O平均鍵長為0.243 nm；Svenum等[22]使用密度泛函理論計算得到C-S-H中Si—O平均鍵長為0.165 nm，Ca—O鍵長為0.23～0.25 nm。由此可見，本文所建模型的Ca—O和Si—O鍵長與實驗數據和密度泛函理論數據一致，說明模型具有可靠性。

不同濕度C-S-H的RDF峰值出現在同一位置，說明濕度并不影響C-S-H中Si—O、Cas—O、Caw—O的鍵長。雖然峰值出現在同一位置，但是峰值大小卻存在差異。隨著濕度的增加，Si—O、Cas—O、Caw—O的RDF峰值不斷增大，說明水分子的增加使得Si、Ca原子在近程范圍內集聚的氧原子數更多。

圖3為三種原子在z方向的原子濃度分布。取4～6 nm的部分晶胞進行分析對比。從圖3可以看出，C-S-H在z方向呈現分層結構，以鈣硅鏈和層間區域交替出現。Si原子和Cas原子組成鈣硅鏈，Caw原子主要分布在兩條鈣硅鏈的層間區域。每個Cas和Caw原子峰值的出現，都伴隨著兩個Si原子峰值的出現，且

圖3 不同濕度C-S-H超胞Si、Cas、Caw原子濃度分布Fig.3 Density profiles of Si， Cas， Caw atoms in C-S-H supercell with different humidity

Cas原子的峰值出現在Si原子的兩個峰值中間，這說明鈣硅鏈由兩條硅鏈組成，且Cas原子分布在兩條硅鏈中間。這與圖1中C-S-H模型所顯示的結構一致。雖然不同濕度的C-S-H的結構一致，但是隨著濕度的增加，層間區域的寬度從0.3 nm增加到0.4 nm再增加到0.47 nm。這說明水分子的進入增大了C-S-H的層間距離，使得C-S-H的分層結構更加明顯。

圖4 層間連接示意圖[6]Fig.4 Schematic diagram of interlayer connection[6]

將圖3揭示的規律，采用圖4表達，圖4為C-S-H的層間連接示意圖[6]。(1)濕度增加會導致C-S-H層間距離變大，體積增加，密度減小，受力時更容易變形。(2)濕度為0%時，硅氧四面體中的O原子與Ca原子跨層連接上下兩個鈣硅層，隨濕度增加，首先層間水分子的O原子與Ca原子成鍵，當濕度進一步增大時，層間又有更多水分子之間的氫鍵，使得相鄰的鈣硅層間較強作用的Si—O鍵和Ca—O鍵逐漸被相對較弱的氫鍵分隔，C-S-H的強度也被削弱。因此，濕度會對C-S-H的力學性能產生重要的影響。

采用1.2節所介紹的方法得到不同濕度的C-S-H的彈性常數,見表1。表中KVoigt、KReuss、KHill分別代表Vogit、Reuss、Hill三種方法所計算的體積模量，GVoigt、GReuss、GHill分別代表Vogit、Reuss、Hill三種方法所計算的剪切模量，Ex、Ey、Ez分別代表三個方向的楊氏模量。從表中可以看到，隨著C-S-H濕度的增加，其體積模量、剪切模量以及三個方向的楊氏模量均不同程度的下降，這說明水分子的侵入會降低C-S-H的彈性性質。Pellenq等[13]通過分子動力學模擬得到濕度為100%時C-S-H的體積模量和剪切模量分別為51 GPa和24 GPa；Abdolhosseini等[23]模擬得到濕度為100%時C-S-H的楊氏模量為67～80 GPa；Hajilar等[24]模擬濕度為100%的C-S-H的體積模量為56 GPa，剪切模量為27 GPa，楊氏模量在69 GPa左右。通過以上對比說明本文所計算的彈性常數具有一定的可靠性。

表1 不同濕度C-S-H超胞的彈性常數Table 1 Elastic constants of C-S-H supercell with different humidity

2.2 不同濕度C-S-H的應力應變關系

圖5為不同濕度的C-S-H單軸拉伸時的應力應變曲線。從圖中可以看出，C-S-H在拉伸時，初始階段應力隨著應變近似線性增加；當應力達到峰值時，應力開始緩降，即出現軟化特性。應力應變關系在兩個方向上存在很大的差異。沿y軸拉伸時，濕度為0%和50%的C-S-H并沒有明顯的軟化特性，即使應變達到0.8時，應力也始終不為零，而濕度為100%的C-S-H在應變達到0.6時，應力減小為0，說明C-S-H已經被拉斷；沿z軸拉伸時，濕度為0%的C-S-H不會被拉斷，而濕度為50%和100%的C-S-H在應變分別達到0.55和0.4時會被拉斷。這說明隨著濕度的增加，C-S-H拉伸時的力學性能和變形能力不斷下降，并且對z方向的影響要大于y方向。沿y方向拉伸時，濕度為0%、50%和100%的C-S-H的彈性極限應變分別為0.048、0.048、0.040，沿z方法拉伸時均為0.040。在彈性階段內進行線性擬合，根據斜率計算y和z方向的彈性模量。三種濕度的C-S-H在y方向的彈性模量分別為77.6 GPa、56.6 GPa、50.5 GPa，z方向的彈性模量分別為70.4 GPa、55.5 GPa、48.0 GPa。顯然隨著濕度的增大，受拉彈性模量不斷減小。在此需要說明的是該方法得到的受拉彈性模量與Materials Studio軟件提供的方法所得計算結果(見表1)存在一定差異，這是由計算原理不同造成的，但濕度影響規律一致。

圖5 C-S-H單軸拉伸應力應變關系曲線Fig.5 Stress-strain relation curves of C-S-H during uniaxial tension

如圖6所示，C-S-H單軸壓縮時的應力應變關系在兩個方向上也存在不同。沿y軸壓縮時，初始階段應力隨應變近似線性增加，三種濕度的極限彈性應變分別為0.048、0.048、0.040。濕度為0%的C-S-H在應力達到峰值后，應力不會減小，反而緩慢增加，表現出強化特性，濕度為50%和100%的C-S-H在應力達到峰值后開始緩慢下降，表現出軟化特性，最后趨于穩定；沿z軸壓縮時，初始階段應力和應變近似呈線性關系，三種濕度的極限彈性應變均為0.048。濕度為0%和50%的C-S-H會表現出強化特性，濕度為100%的C-S-H會表現出軟化特性，最后趨于穩定。說明濕度的增加會降低C-S-H壓縮時的力學性能，且濕度對抗壓能力的影響要小于對抗拉能力的影響。三種濕度的C-S-H在y方向的受壓彈性模量分別為88.9 GPa、82.3 GPa 80.7 GPa，在z方向的受壓彈性模量分別為57.6 GPa、55.3 GPa、49.5 GPa，受壓彈性模量隨濕度的增加而降低。其規律與Materials Studio軟件提供方法計算結果(見表1)相同，但值存在差異。另外，由于受力變形特性不同，與拉伸時結果也存在差異，說明C-S-H拉、壓彈性模量是不同的。

圖6 C-S-H單軸壓縮應力應變關系曲線Fig.6 Stress-strain relation curves of C-S-H during uniaxial compression

圖7為不同濕度C-S-H剪切時的應力應變關系曲線。C-S-H在受剪時，初始階段應力與應變近似呈線性關系，三種濕度的極限彈性應變均為0.035。當剪應力達到峰值后，應力呈鋸齒狀波動并緩慢增加，表現出強化特性。另外，隨著濕度的增加，剪應力在不斷下降。說明濕度的增加會降低C-S-H剪切時的力學性能，但影響程度也小于對抗拉性能的影響。根據彈性階段的應力應變關系得到三種濕度的C-S-H剪切模量分別為19.7 GPa、17.4 GPa、13.5 GPa，表明隨著濕度的增加，剪切模量逐漸降低。其規律與Materials Studio軟件提供方法計算結果(見表1)相同，但值存在差異。

圖7 C-S-H剪切應力應變關系曲線Fig.7 Stress-strain relation curves of C-S-H under shearing

以上通過對應力應變曲線的簡單分析可以看出，濕度對C-S-H力學性能存在影響，同時，C-S-H在不同的受力狀態下表現出不同的力學性能和變形性能。為了進一步分析濕度對不同受力狀態下C-S-H力學性能的影響程度，表2列出了不同濕度的C-S-H在三種受力狀態下的應力應變特征參數。從峰值應力的降低程度來看，當C-S-H的濕度從0%增加到50%和100%，在受拉狀態下，y方向分別降低了4.90%和9.92%，z方向分別降低了6.56%和26.65%；在受壓狀態下，y方向分別降低了2.88%和4.01%，z方向分別降低了1.79%和9.59%；受剪時，分別降低了7.34%和21.96%。說明濕度對C-S-H抗拉與抗剪強度的影響較大，對抗壓強度的影響較小。

應力應變曲線下的面積表示C-S-H的應變能密度，反映了C-S-H破壞前吸收的能量大小和變形性能。從應變能密度的降低程度來看，當C-S-H的濕度從0%增加到50%和100%，在受拉狀態下，y方向分別降低了11.18%和53.11%，z方向分別降低了42.88%和68.29%；在受壓狀態下，y方向分別降低了11.06%和24.26%，z方向分別降低了10.55%和20.66%；受剪時，分別降低了20.16%和21.49%。說明濕度對C-S-H拉伸時的變形性能影響最大，其次是剪切，對壓縮時變形性能的影響程度最小。

表2 不同受力狀態下C-S-H的應力應變特征參數Table 2 Characteristic parameters of stress-strain of C-S-H under different stress states

3 結 論

(1)濕度增加使得C-S-H中Si、Ca原子近程范圍內集聚的氧原子數更多，同時層間距離不斷增大，分層更加明顯；另外濕度會降低C-S-H的體積模量、剪切模量以及楊氏模量。

(2)濕度的增加會降低C-S-H單軸拉伸、單軸壓縮和剪切時的力學性能和變形性能，對沿鏈狀方向(y向)力學性能的影響弱于垂直鏈狀方向(z向)。

(3)濕度對抗拉與抗剪強度的影響較大，對抗壓強度的影響較??；對拉伸時的變形性能影響最大，其次是剪切，對壓縮時變形性能的影響最小。