SiO2氣凝膠分子動力學模擬研究進展

楊 云，史新月，吳紅亞，秦勝建，張光磊

(石家莊鐵道大學材料科學與工程學院，石家莊 050043)

0 引 言

氣凝膠材料因其三維納米骨架結構，具有低密度、高孔隙率、高比表面積的特性，符合目前追求輕質、堅固、經濟的超輕工程材料的發展趨勢，近年來得到廣泛關注。制備氣凝膠的材料選擇上有很大的靈活性，Kistler等[1]早在二十世紀三十年代就通過超臨界干燥法成功合成了二氧化硅、氧化鋁和氧化鎢氣凝膠。二氧化硅氣凝膠是目前研究最成熟，使用最廣泛的無機氧化物氣凝膠。其最早的化學前驅體是偏硅酸鈉[2]，目前逐步被硅醇氧化物取代，例如正硅酸甲酯、正硅酸乙酯和甲基三甲氧基硅烷。這些單體縮聚形成多孔凝膠，再經過超臨界干燥得到最終的氣凝膠結構。

由于獨特的結構和性能，二氧化硅氣凝膠在工程方面應用較為廣泛。二氧化硅氣凝膠可以識別、捕捉粒子，被用于研制切倫科夫輻射器[3]、被動式太陽能收集設備[4]以及捕捉太空塵埃[5]，甚至被用于大量納米粒子的捕獲、存儲和釋放[6-7]。其大規模、低成本的生產在隔熱領域也有很大發展。Wittwer等將二氧化硅氣凝膠用于被動式太陽能窗[8]、窗戶隔熱[9]以及航天運載工具上的低溫隔熱系統[10]。

二氧化硅氣凝膠在眾多領域都顯示出廣闊的應用前景，但是由于氣凝膠中無機網格骨架的脆性以及高溫穩定性差等原因，其實際應用受到限制。目前，從原子結構角度探究氣凝膠的組成成分、分形結構特征及熱力學性能之間的關聯，分析得到熱力學性能更優的氣凝膠結構并指導實驗，已經逐漸成為國內外的研究熱點。分子動力學方法(MD)可以直觀地呈現二氧化硅氣凝膠的內部成鍵方式和分形結構特征，并且能夠模擬二氧化硅氣凝膠內部結構與力學性能的關系。Ng[11]和Yeo[12]等應用MD預測了氧化硅氣凝膠中密度與熱傳導系數的關系，模擬結果與實驗相吻合；Rivas[13]和Lei[14]等應用MD模擬了氣凝膠密度與彈性模量之間的關系；Ferreiro等[15]應用MD模擬了二氧化硅氣凝膠中塑性變形與斷裂過程中的微結構演化規律；Patil等[16]應用MD揭示了氧化硅氣凝膠、石墨烯增強的二氧化硅氣凝膠和納米玻璃纖維增強的二氧化硅氣凝膠的力學行為與微結構演變規律。接下來本文將從二氧化硅氣凝膠的結構特性、勢函數選取、模型建立以及性能等方面進行介紹。

1 二氧化硅氣凝膠的特性

在標準狀況下空氣的密度大約是0.001 2 g/cm3，而二氧化硅氣凝膠由于其多孔性，其密度可以低至0.003 g/cm3，僅僅是空氣密度的三倍左右[17]。

表征二氧化硅氣凝膠力學性能的方法有很多種。Cross等[18]通過凝膠法測量聲速，得到在密度為0.05～2 g/cm3的氣凝膠中，其彈性模量和剪切模量的密度冪律指數分別為3.17±0.11和3.25±0.15。只有在中等密度0.158 g/cm3時，聲速法和單軸壓縮試驗的方法測量的楊氏模量才會表現出一致性。進一步的單軸壓縮測試發現，彈性模量隨密度的變化因不同的處理條件而不同。在密度范圍為0.1～0.6 g/cm3之間，酸催化的氣凝膠的冪律指數為3.37±0.13，而堿催化的氣凝膠為2.85±0.17，預水解的氣凝膠為3.77±0.21[19]。對密度范圍為0.1～0.4 g/cm3之間的氣凝膠進行梁彎曲試驗，確定彈性模量隨密度的變化，得出其冪律指數為3.8[20]。二氧化硅氣凝膠的低密度和高孔隙率，在外力作用下容易發生脆性斷裂。大量的研究表明，改變前驅體的種類可以獲得較好的韌性。

雖然多孔性不利于氣凝膠的力學性能，但正是由于其多孔結構，二氧化硅氣凝膠具有異常低的熱導率。氣凝膠材料中的熱量傳遞主要有三條途徑：固態傳熱、氣態傳熱以及輻射傳熱[21]。降低熱導率的因素有很多，例如忽略不同傳熱模式之間的耦合；氣孔內氣體自由運動受束縛而限制了對流傳熱；氣凝膠骨架內的局部激勵受到有效限制等[22]。特別是Fricke等[23]發現了固體熱導率(λs)隨密度(ρ)的冪律變化關系：

λs=Cρα

(1)

式中：C隨氣凝膠不同而變化；在密度為0.3～1 g/cm3之間α為1.6，在密度為0.07～0.3 g/cm3之間α為1.2。

2 MD模擬下二氧化硅氣凝膠模型的建立

2.1 勢函數

在MD模擬中，原子間電位起著至關重要的作用，包括成鍵方式和與鄰近原子間的相互作用。常用的勢函數有BKS勢、Vashishta勢和Tersoff勢等。

BKS勢(van Beest, Kramer and van Santen)[24-25]由Carré等提出，其數學公式為：

(2)

式中，第一部分和第二部分分別表示原子的短程引力和斥力，最后一部分指長程庫侖相互作用力。rij表示i和j原子間距，qi和qj表示原子電荷，Aij、Bij和Cij表示不同對原子i和j的常數。BKS勢可以準確地描述各種硅晶型的結構和力學性能[26]，其缺陷在于高淬火溫度下，原子之間會過于接近。這個缺陷可以通過加入“24-6” Lennard-Jones勢[27]來克服，其數學公式如下：

(3)

式中，參數εij和σij取值如下：εSi-Si=13.20 eV，εSi-O=1.12×10-2eV，εO-O=4.78×10-4eV；σSi-Si=0.04 nm，σSi-O=0.135 nm，σO-O=0.22 nm[12]。

BKS勢包含一個長程庫侖相互作用項，因此計算大系統的力非常復雜。Carré等引入了截止庫侖長程相互作用的wolf方法[28]，如下：

(4)

(5)

(6)

其中：

(7)

(8)

(9)

式中，α，β分別表示硅原子和氧原子，e表示電子，參數rc,sh=0.55 nm，γw=γsh=0.05 nm，參數qα，qβ，Aαβ，Bαβ，Cαβ取值可參考文獻[24]。Vw(r)表示施加一個有限的截止，在長程庫侖相互作用下rc,w=1.217 nm，從而避免了力在r=rc,w處不連續。

除了BKS勢，一些研究應用Vashishta勢[29]來模擬二氧化硅氣凝膠的結構和性能[30-31]。Vashishta勢的總能量表示如下：

U=U2(rij)+U3(rij,rjk,rik)

(10)

兩體相互作用項為：

(11)

式中：Hij和ηij分別表示空間斥力的強度和指數；Zi和Zj代表有效電荷；αi代表電子極化率；r4s是偶極相互作用的衰變長度。

三體相互作用項為：

(12)

式中，Bijk是相互作用強度。鍵的拉伸效應由以下函數表示：

(13)

式中r0是相互作用的截斷距離。鍵的彎曲效應由以下函數表示：

(14)

式中θijk是以原子i為頂點的原子i、j、k的角度。

為了更準確地模擬每個原子所處的化學環境，需要解釋多體相互作用的關系。Tersoff勢[32]可以通過改變每個鍵的強度來解決這一問題，并將二氧化硅氣凝膠參數化，更好地預測其熱性能和力學性能[12,14,33]。其數學公式如下：

U=fc(rij)[VR(rij)+bijVA(rij)]

(15)

式中：VR和VA分別表示排斥項和吸引項；fc是限制短程相互作用的力；bij是根據i和j原子的配位而單調遞減的鍵序;排斥項和吸引項的表達式如下：

VR(rij)=Aijexp(-γijrij)

(16)

VA(rij)=Bijexp(-μijrij)

(17)

(18)

截斷函數如下：

(19)

(20)

式中，R和S表示相互作用極限。鍵序函數項如下：

(21)

(22)

(23)

式中，χij表示異極鍵的增強和減弱，Ai，Aj，Bi，Bj等表示勢函數中的參數，其取值可參考文獻[12]。

2.2 二氧化硅氣凝膠模型的建立

LAMMPS[34]是進行MD的開源仿真軟件，由美國Sandia國家實驗室開發，可以免費獲取使用，使用者可以根據需要自行修改源代碼。LAMMPS可以支持包括氣態、液態和固態相形態下、各種系綜、百萬級的原子分子體系，并支持多種勢函數，具有良好的并行擴展性。

二氧化硅氣凝膠的分形結構由Kieffer方法[35]產生。生成二氧化硅氣凝膠的步驟為：(1)建立模擬盒子，以β-方石英晶體建立含有a×b×c個單胞的SiO2晶體盒子；(2)模型采取三維周期性邊界條件，運用共軛梯度法進行能量最小化，得到在絕對零度下穩定的晶體結構；(3)在NPT(恒定壓力和溫度)系綜下，高溫(一般為5 000 K或6 000K)平衡一段時間，將晶體熔化，然后淬火至300 K，得到非晶態氧化硅；(4)在NVE(恒定體積和能量)系綜下，溫度為300 K，進行結構弛豫，使拉伸-松弛階段不斷重復，在松弛階段通過Si-O鍵的逐漸斷裂形成孔隙，得到一定密度的氣凝膠；(5)將氣凝膠平衡一段時間，使其結構在恒定的壓力和體積下處于穩定狀態。以上步驟可以根據需要進行改動。圖1為二氧化硅氣凝膠模型各階段的示意圖，其中，(a)為晶體SiO2，(b)為非晶二氧化硅，(c)為二氧化硅氣凝膠。

圖1 二氧化硅氣凝膠模型建立過程各階段的示意圖Fig.1 Schematic diagram of each stage of silica aerogel modeling process

負壓破裂法是另一種常用的二氧化硅氣凝膠建模方法。Kieffer和Angell通過MD對二氧化硅氣凝膠的滲透結構進行建模[34]。在“負壓破裂”方案中，等速分離的非晶二氧化硅樣品立即逐步膨脹，直到達到滲透閾值。所有的鍵在每一個擴展步驟中都被同時拉伸到超出其平衡10%～35%的長度，隨后結構被松弛以消除過多的累積應力。Nakano等用負壓破裂的方法和Vashishta勢建立了包含40 000 個原子的二氧化硅氣凝膠模型[30]。Ng等應用負壓破裂方法和BKS勢建立了3 000個原子的二氧化硅氣凝膠模型[11]。

除以上兩種方法外，一種改進的負壓破裂法也被應用到氣凝膠結構建模中[13]，包括膨脹、加熱和淬火[12,36]。(1)膨脹：通過擴大體積的方法使β-方石英晶體晶格中原子膨脹到所需的密度。(2)加熱：將整個系統加熱到3 000 K來擴散和聚集。(3)淬火：將溫度降低到2 500 K和0 K之間來松弛。這種方法在結構特性上與試驗結果有很好的一致性，例如，分形維數、徑向分布函數和模擬散射的結果。

3 MD模擬下二氧化硅氣凝膠的表征

3.1 結構表征

氣凝膠模型具有特殊結構和特性，如孔徑大小分布、分形維數、表面積、鍵長和鍵角分布[37]。徑向分布函數g(r)可以用來區分晶體和非晶體，表征二氧化硅氣凝膠的徑向分布函數可以得到Si-O鍵、O-O鍵和Si-Si鍵的鍵長。Rivas等[13]通過徑向分布函數圖，得到二氧化硅氣凝膠中Si-O鍵、O-O鍵和Si-Si鍵的鍵長為(0.163±0.003) nm、(0.265±0.003) nm和(0.308±0.003) nm，與實驗上二氧化硅的特征相似。

分形維數df是描述二氧化硅氣凝膠分形特征的重要參數。Kieffer等[35]提出基于總徑向分布函數g(r)的長范圍衰減的方法來計算分形維數，并給出了如圖2所示的徑向分布函數圖，分形維數df表示如下：

(24)

式中，r為原子間的距離。二氧化硅氣凝膠的分形維數根據其加工條件的不同而變化，基本加工條件下的分形維數約為1.8[38]，酸性和中性條件下[39]的分形維數約為2.2～2.4。對比Gon?alves[40]、Kieffer[35]和Rivas[13]的模擬結果，分形維數隨密度的增大而增加，并且接近于線性增長趨勢。Woignier等[41]的實驗分析得出，在溶膠-凝膠過程中，二氧化硅氣凝膠的分形維數根據溶劑酸度的不同在1.8～2.4之間變化，這表明分形維數非常依賴于加工條件。

圖2 不同密度二氧化硅氣凝膠徑向分布函數[35] (插入的對數是用來進行分形維數計算的區域)Fig.2 Radial distribution function of silica aerogel with different density. The insert indicates the zone used to compute the fractal imension[35]

圖3 密度為0.45 g/cm3模型的散射強度與波數的對數圖[17]Fig.3 Logarithmic diagram of scattering intensity and wave number for a model with a density of 0.45 g/cm3[17]

另一種獲得分形維數的方法是通過模擬散射實驗[42-43]。為了進行模擬散射實驗，需要計算散射強度I，對應于樣品所受的不同波長的輻射，用它們的波數q表示。當系統中所有粒子的位置都知道時，計算I的表達式如下[44-45]：

(25)

式中，I是散射強度，I0是強度的參考值，q是波數，rij是兩個原子i和j之間的距離。對于分形結構，散射強度可以通過冪律與波數相關得到，其中指數是分形維數的負值，I∝q-df。Rivas等[13]繪制了I(q)，如圖3所示。

3.2 力學性能

3.2.1 彈性模量

利用MD模擬來研究二氧化硅氣凝膠的力學性能，可以從原子尺度揭示脆性的本質并提出可能增韌的方法。從目前的研究來看，二氧化硅氣凝膠的彈性模量服從冪律縮放，這是分形結構的典型特征。Gon?alves等[40]使用Berendsen恒壓器[46]和Langevin恒溫器[47]對體系進行力學測試，使用加載-卸載循環的應力模式，得到彈性模量隨密度變化的冪律指數為3.84±0.22，與Woignier等的實驗結果吻合較好[20,48]。Campbell等使用Vashishta勢生成密度在1.67～2.2 g/cm3之間的二氧化硅氣凝膠[49]，楊氏模量的冪律指數為3.5±0.2。Berendsen等[46]采用Vashishta勢構建密度在0.23～2.2 g/cm3之間的二氧化硅氣凝膠，并研究其力學性能。楊氏模量和強度的冪律變化指數分別為3.11±0.21和2.53±0.15，與實驗數據吻合良好。Patil等用同樣的方法構建了密度范圍為0.28～2.2 g/cm3的二氧化硅氣凝膠模型[31]，其楊氏模量的冪律指數約為3.25±0.10。表1表明二氧化硅氣凝膠冪律指數與模型大小的相關性，模型越大，冪律指數越大。

表1 模型大小與冪律指數相關性[40]Table 1 Sample size is correlated with power law exponents[40]

3.2.2 抗拉強度

二氧化硅氣凝膠的抗拉強度是表征力學性能的重要參數。Goncalves等[40]對9個不同體積(203～1003nm3)的二氧化硅氣凝膠樣品進行單軸拉伸試驗，結果表明：(1)二氧化硅氣凝膠的抗拉強度是一個與體積有關的性能指標；(2)當體積>403nm3時，抗拉強度與體積成冪律關系，并隨體積的增大而減小，如圖4所示。

圖4 抗拉強度與樣品體積的對數關系[40]Fig.4 Logarithmic relationship between tensile strength and sample volume[40]

圖5 (a)斷裂韌性KIC和(b)應變能釋放率G與裂紋長高比a/h的關系[51]Fig.5 Influence of the crack length to height ratio a/h on the (a) fracture toughness KIC and (b)strain energy release rate G[51]

3.2.3 斷裂韌性

Patil等[51]以507萬個原子為模型，利用Vashishta勢[34]對二氧化硅氣凝膠進行建模，研究了裂紋長高比對斷裂韌性和應變能釋放率的影響。結果如圖5所示。在裂紋擴展初期，斷裂韌性KIC和應變能釋放率G隨裂紋長高比的增大而增大。當比率大于0.2～0.25時，斷裂韌性KIC在所有密度下保持恒定。KIC和G對密度的冪律關系分別為2.02±0.05和1.16±0.04。對于所有密度的氣凝膠，隨著裂紋長高比的增大，斷裂強度減小。以上結果為氣凝膠的多尺度模擬奠定了基礎。

3.3 傳熱性能

3.3.1 熱導率計算

氣凝膠材料中的熱量傳遞主要有三條途徑：固態傳熱、氣態傳熱以及輻射傳熱。氣凝膠材料中總的有效熱導率可表示為[52]：

λtotal=λs(T)+λg(T,pg)+λr(T)

(26)

式中：λtotal、λs(T)、λg(T,Pg)、λr(T)分別表示氣凝膠的總熱導率、固相熱導率、氣相熱導率和輻射熱導率。因為氣凝膠的孔隙率可以達到90%以上，這種特殊結構決定了其具備非常好的隔熱性能。

固相熱傳導主要取決于氣凝膠內部的骨架結構，其傳熱過程可看成是聲子在氣凝膠骨架中的傳播，可通過聲子擴散的德拜模型進行描述：

(27)

式中：λ0、ρ0、CV、lph、V0分別表示氣凝膠骨架熱導率、骨架密度、骨架定容比熱容、聲子平均自由程以及波的傳播速度。

氣凝膠材料中氣相傳熱主要包括孔結構的氣相熱傳導和氣相熱對流，由于氣凝膠材料的孔徑大部分處于介孔范圍內，小于空氣分子平均自由程，因此氣體分子在氣凝膠孔內傳遞時會與孔骨架相撞。氣凝膠中氣相傳熱可用Kaganer等提出的公式來計算[53]：

(28)

式中：π、λp0、β、lg、D、Kn分別表示氣凝膠孔隙率、無對流自由氣體的熱導率與氣體種類和溫度相關的系數、氣體分子平均自由程、孔徑、Knudsen數。

輻射傳熱是一種非接觸式的傳熱方式，不需要任何介質，可以看成是光子與氣凝膠材料的相互作用。在溫度較低時，輻射傳熱不明顯，可忽略不計，但是當外界溫度高于600 ℃時，輻射傳熱在總的熱導率中將占主導。輻射傳熱可通過下式進行計算[54]：

(29)

圖6 熱導率與密度的關系[11]Fig.6 Relation between thermal conductivity and density[11]

式中：n、Tr、ρ、e分別表示氣凝膠有效折光系數、平均溫度、密度以及消光系數。

3.3.2 傳熱模擬

MD模擬也被廣泛用于二氧化硅氣凝膠固體骨架中納米尺度熱傳輸的表征。一系列的MD模擬使用非平衡分子動力學(NEMD)的方法研究了二氧化硅氣凝膠的熱導率[55]。Ng等[11]應用BKS勢建立二氧化硅氣凝膠模型，并用反向非平衡動力學的方法(RNEMD)測定其熱導率[56]。在648個原子的體系中，熱導率為(1.20±0.07) W/(m·K)，與之前MD結果一致[57-58]；3 000個原子的系統熱導率為1.39 W/(m·K)，接近實驗值1.40 W/(m·K)[12]。Ng等應用負壓破裂的方法得到密度在0.32～0.99 g/cm3之間的氣凝膠,并使用RNEMD測定熱導率，如圖6所示。結果表明，熱導率的冪律指數為1.01，與實驗值有一些差距。

Yeo等[12]應用重新參數化的Tersoff勢，結合膨脹、加熱和淬火方案建立二氧化硅氣凝膠模型。結果表明，熱導率隨密度的冪律指數為1.61，與實驗值保持一致。該研究進一步表明，二氧化硅氣凝膠顆粒之間的熱阻[50]對熱導率的測定有一定的影響。

4 結 論

二氧化硅氣凝膠的低密度、低熱導率和高孔隙率等優異性能使其在航空航天、交通運輸和建筑保溫等領域有極為廣闊的應用前景，但是其低強度和高脆性的特點限制了其大規模應用。了解其內部結構與力學性能的關聯機制是設計二氧化硅氣凝膠結構并改善其力學性能的關鍵。

分子動力學方法可以從原子尺度分析二氧化硅氣凝膠的微觀結構并預測其熱力學性能。利用分子動力學方法，選取合適的勢函數可以重現二氧化硅氣凝膠分形結構模型。通過統計徑向分布函數和分形維數可以得到二氧化硅氣凝膠的鍵長分布和分形結構特征。彈性模量、抗拉強度和熱導率是衡量二氧化硅氣凝膠熱力學性能的重要參數。應用分子動力學方法可以模擬得到二氧化硅氣凝膠彈性模量、抗拉強度和熱導率與密度之間的冪律關系。本項工作有助于揭示二氧化硅氣凝膠的熱力學性能、絕熱能力以及微結構演變規律，對于提升二氧化硅氣凝膠熱力學性能并擴大其應用領域有重要的指導意義。