高壓下HBT 晶體的彈性性質(zhì)

李 佐，張鳳玲，廖大麟

（貴州工程應(yīng)用技術(shù)學(xué)院理學(xué)院，貴州 畢節(jié) 551700）

含能材料的彈性性質(zhì)包括彈性常數(shù)、體彈模量、楊氏模量和剪切模量等。它反映材料在常溫、靜荷載作用下的宏觀力學(xué)性能，不但決定材料對施加應(yīng)力的響應(yīng)方式，還能反映含能材料的穩(wěn)定性等[1]。目前，含能材料的彈性性質(zhì)得到了人們的廣泛關(guān)注，并對RDX[2]、TATB[3]、HMX[4]、FOX-7[5]、CL-20[6]等含能材料開展了研究。5, 5′-雙四唑肼(5, 5′-Hydrazinebistetrazole，HBT)晶體是一種著名的單質(zhì)富氮、高能、鈍感含能材料[7]。HBT 含有83.7%的氮元素，爆炸燃燒后釋放的氣體主要為氮?dú)?、二氧化碳和水，因此被稱為“綠色”含能材料[8-11]。本研究采用基于第一性原理的準(zhǔn)諧振近似方法，計(jì)算固態(tài)HBT 晶體在0～120 GPa 高壓下的彈性性質(zhì)，同時(shí)分析其彈性的各向異性特征。

1 理論與計(jì)算方法

1.1 HBT 晶體結(jié)構(gòu)優(yōu)化方法

采用基于第一性原理的Quantum ESPRESSO 軟件包中的PWSCF 模塊[12]計(jì)算電子結(jié)構(gòu)的總能量。廣義梯度近似（GGA）下交換相關(guān)勢采用Perdew-Burke-Ernzerhof （PBE）形式[13-14]，原子贗勢平面波基組采用超軟贗勢加核修正的C.pbe_v1.2.uspp.F.UPF，N.pbe_v1.2.uspp.F.UPF，H.pbe_v1.2.uspp.F.UPF 形式[15]，幾何優(yōu)化采用BFGS 算法[16]。為增加計(jì)算精度，平面波基函數(shù)的截?cái)嗄苋?44 eV，布里淵空間采用Monkhorst-Pack 方法[17]，k 點(diǎn)方案為6 × 6 × 3，總能自洽，離子位移最大值和應(yīng)力最大值都能達(dá)到設(shè)置的收斂精度。在高壓下幾何結(jié)構(gòu)優(yōu)化的基礎(chǔ)上，利用準(zhǔn)諧振近似方法計(jì)算HBT 晶體的彈性性質(zhì)，詳細(xì)的計(jì)算過程在軟件Thermo_pw 中體現(xiàn)[18]。

1.2 彈性常數(shù)

固態(tài)物質(zhì)彈性理論計(jì)算中引入應(yīng)變εj（j = 1, 2, ···, 6），晶體的一個(gè)應(yīng)力用σi表示。當(dāng)應(yīng)力較小時(shí)，σi和εj滿足[19]

式中：系數(shù)Cij為彈性常數(shù)（或者彈性勁度），矩陣S 為矩陣C 的逆矩陣，Sij為彈性柔順系數(shù)。初始構(gòu)型采用Klap?tke 等[7]獲得的HBT 晶體的實(shí)驗(yàn)結(jié)構(gòu)，它屬于單斜結(jié)構(gòu)的分子晶體，點(diǎn)群為C2/c 空間群，其原胞晶系屬于三斜晶系，對稱性最低，包含最多的矩陣元。彈性常數(shù)矩陣C 是一個(gè)對稱矩陣，共有21 個(gè)獨(dú)立矩陣元。計(jì)算中采用原胞結(jié)構(gòu)，彈性常數(shù)矩陣表達(dá)式為

考慮到二階彈性常數(shù)與能量滿足如下關(guān)系[20]

因此，只要知道穩(wěn)定構(gòu)型下的晶體體積V0和能量Ec，就可以獲得彈性常數(shù)。

2 結(jié)果與討論

2.1 晶體結(jié)構(gòu)參數(shù)

基于密度泛函理論的第一性原理，得到零溫零壓下HBT 晶體的晶格參數(shù)，如表1 所示，晶體幾何結(jié)構(gòu)如圖1 所示。從表1 可以看出，晶格參數(shù)b、c 和V 的計(jì)算值均高于實(shí)驗(yàn)結(jié)果，而a 和β 卻低于實(shí)驗(yàn)結(jié)果。在目前的第一性原理方法中，晶格參數(shù)的計(jì)算精度取決于原子贗勢和計(jì)算方法的選擇，因此，無法做到與實(shí)驗(yàn)結(jié)果完全符合。此外，采用5 GPa 間隔優(yōu)化HBT 晶體在0～120 GPa 壓強(qiáng)范圍內(nèi)的晶體結(jié)構(gòu)。優(yōu)化后獲得了穩(wěn)定的晶體構(gòu)型，壓強(qiáng)-體積比（p-V/V0）曲線與實(shí)驗(yàn)結(jié)果的比較如圖2 所示。從圖2 可以看出，在0～26 GPa 壓強(qiáng)范圍內(nèi)，理論計(jì)算結(jié)果與Ciezak-Jenkins 等[21]的實(shí)驗(yàn)結(jié)果具有相同的變化趨勢，并且在1 和22 GPa 處曲線相交，說明本理論計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果有一定的可比性。同時(shí)，給出了HBT 晶體在壓強(qiáng)高于25 GPa 時(shí)的體積變化情況。理論計(jì)算顯示，隨著壓強(qiáng)的增大，晶體體積變化受壓強(qiáng)的影響越來越明顯，60 GPa以上時(shí)壓強(qiáng)與體積比呈線性關(guān)系。

表 1 HBT 晶體晶格參數(shù)的計(jì)算值和實(shí)驗(yàn)值Table 1 Calculated lattice parameters of HBT crystal along with experimental data

圖 1 優(yōu)化后的HBT 晶體幾何結(jié)構(gòu)Fig. 1 Geometric structure of optimized HBT crystal

圖 2 理論和實(shí)驗(yàn)得到的壓強(qiáng)-體積比關(guān)系Fig. 2 Pressure as a function of volume ratio V/V0 from theoretical and experimental data

2.2 高壓彈性

通過高壓下結(jié)構(gòu)優(yōu)化，得到不同壓力下穩(wěn)定的HBT 晶體結(jié)構(gòu)，在此基礎(chǔ)上計(jì)算彈性系數(shù)和彈性模量。三斜結(jié)構(gòu)的彈性常數(shù)矩陣包含21 個(gè)獨(dú)立矩陣元。Born[22]d 給出了晶體材料力學(xué)穩(wěn)定性的充分必要條件為彈性常數(shù)矩陣C 正定。按照線性代數(shù)知識可知，對稱矩陣C 正定的充要條件為C 的各階順序主子式都為正，即

式(5)～式(8)給出了三斜晶體穩(wěn)定性的充要條件。四階以上行列式的表達(dá)式比較復(fù)雜，一般采用直接帶入行列式元素來計(jì)算驗(yàn)證。圖3(a)、 圖3(b)、圖3(c)、圖3(d)給出了21 個(gè)彈性常數(shù)隨壓強(qiáng)的變化關(guān)系。從圖3 可以看出，C15、C56、C24、C46隨壓強(qiáng)的變化比較劇烈，其中C15、C56和C24先減小后增大，而C46則先增大后減小隨后又增大。將零溫零壓下晶體的彈性常數(shù)代入式(5)～式(8)，驗(yàn)證后發(fā)現(xiàn)符合穩(wěn)定性條件。對高壓下晶體進(jìn)行驗(yàn)證，發(fā)現(xiàn)也符合穩(wěn)定性條件，說明高壓下HBT 晶體是穩(wěn)定的。

計(jì)算體彈模量、剪切模量的普適公式為[23]

根據(jù)Voigt-Reuss-Hill 均值方法，有

再根據(jù)B、G 可以獲得體彈模量E 和泊松比μ

依據(jù)式(9)～式(11)，可以獲得體彈模量、剪切模量和楊氏模量隨壓強(qiáng)的變化關(guān)系，如圖4 所示。從圖4 可以看出，隨著壓強(qiáng)的增大，3 種模量的數(shù)值也增大，其中體彈模量在15 GPa 以上時(shí)表現(xiàn)出與壓強(qiáng)的線性關(guān)系。而3 種模量中，楊氏模量的變化最顯著，高壓下的數(shù)值也最大。此外，根據(jù)Pugh[24]的理論，B/G 可作為鑒別晶體韌脆性的標(biāo)準(zhǔn)，即當(dāng)B/G＜1.74 時(shí)，晶體材料表現(xiàn)為脆性；反之，晶體材料表現(xiàn)為韌性。因此，根據(jù)式(11)中泊松比 μ與B/G 的關(guān)系來計(jì)算B/G 的數(shù)值，并作出B/G 與壓強(qiáng)的關(guān)系曲線，如圖5 所示。隨著壓強(qiáng)增大，B/G 逐漸增大，與圖4 的變化趨勢一致。當(dāng)壓強(qiáng)低于85 GPa 時(shí)，B/G＜1.74，材料表現(xiàn)出脆性；而當(dāng)壓強(qiáng)高于85 GPa 時(shí)，HBT 晶體表現(xiàn)出韌性。類似的層狀分子晶體TATB 在0～50 GPa 時(shí)表現(xiàn)為韌性[25]，而同為層狀分子晶體的HBT 在0～85 GPa 下卻表現(xiàn)為脆性，說明HBT 晶體對壓力的響應(yīng)較為敏感。

圖 3 彈性常數(shù)隨壓強(qiáng)的變化Fig. 3 Variations of elastic constant with pressure

圖 4 體彈、剪切和楊氏模量隨壓強(qiáng)的變化關(guān)系Fig. 4 Variations of bulk modulus B, shear modulus G and Young’s modulus E with pressure

圖 5 體彈模量和剪切模量之比與壓強(qiáng)的關(guān)系Fig. 5 Pressure dependence of the ratio of bulk modulus and shear modulus

同時(shí)，根據(jù)以下公式可以得到HBT 晶體材料的壓縮波聲速 vp、剪切波聲速 vs及 平均聲速vm

式中：h 為普朗克常數(shù)，kB為玻爾茲曼常數(shù)，n 為原子數(shù)，NA為阿伏加德羅常數(shù)，M 為相對分子質(zhì)量。圖6 給出了壓縮波聲速 vp、剪切波聲速 vs、平均聲速 vm和德拜溫度 Θ隨壓強(qiáng)的變化關(guān)系?？梢姡曀匐S壓強(qiáng)的增大而增大，其中壓縮波聲速 vp最大，且隨壓強(qiáng)變化最明顯。在0 GPa 下得到的德拜溫度為975 K，而120 GPa 時(shí)則為1 820 K，壓強(qiáng)對德拜溫度的影響明顯。

圖 6 不同壓強(qiáng)下HBT 晶體的聲速和德拜溫度Fig. 6 Compressional wave velocity, shear wave velocity and averaged wave velocity of HBT crystal under different pressures

2.3 高壓彈性各向異性

彈性性質(zhì)的各向異性是材料力學(xué)性能通常受微裂紋和晶格畸變影響的重要條件，因此，研究彈性性質(zhì)的各向異性對于提高材料的力學(xué)性能有著重要的意義[27]。目前有3 種處理方案。第1 種方案是Chung 等[28]提出的體彈模量、楊氏模量和剪切模量分?jǐn)?shù)比各向異性的概念，定義為

式(16)～式(18)分別用于表征材料的體彈模量、楊氏模量和剪切模量的各向異性程度。據(jù)此可以計(jì)算出模量各向異性分?jǐn)?shù)比與壓強(qiáng)的關(guān)系，如圖7 所示。對于彈性各向同性材料，AB、AE和AG都為零，而AB、AG為1 表示最大可能的彈性各向異性。從圖7 可以看出，AB、AE和AG均大于零而小于1，說明在0～120 GPa 的壓強(qiáng)范圍內(nèi)HBT 晶體的彈性模量表現(xiàn)出各向異性。此外，體彈模量分?jǐn)?shù)比只有很小的波動，楊氏模量和剪切模量則表現(xiàn)出相同的變化趨勢，即先增大后逐漸減小，說明HBT 晶體在高壓下的各向異性程度減弱。

圖 7 模量分?jǐn)?shù)比與壓強(qiáng)的關(guān)系Fig. 7 Pressure dependence of the fractional ratio of modulus

第2 種方案是Ranganathan 等[29]提出的全局彈性各向異性系數(shù)（AU）的概念，用來描述材料彈性模量的各向異性，具體表達(dá)式如下

圖8 的縱軸表示0～120 GPa 壓強(qiáng)范圍內(nèi)的全局彈性各向異性系數(shù)AU的數(shù)值，其中，最高點(diǎn)表示0 GPa 下的數(shù)值。從結(jié)果看都大于零，因此，HBT 晶體在高壓下表現(xiàn)出各向異性。

此外，式(20)的AC與式(18)的AG是相等的，都是由Chung 等于1967 年提出，均用于表征剪切模量的各向異性。為了比較兩種方案的關(guān)系，作AU-AC曲線?？梢姡瑑煞N方法都很好地描述了HBT 晶體在0～120 GPa 壓強(qiáng)范圍彈性模量的各向異性特征。隨著壓強(qiáng)的增大，AU、AC和AU/AC的數(shù)值都在減少，說明彈性模量的各向異性程度在減弱，兩種方案得到的結(jié)論一致。

第3 種方案是Toher 總結(jié)和發(fā)展的Ranganathan 關(guān)于材料彈性各向異性的理論表述[30]

AL隨壓強(qiáng)的變化關(guān)系如圖9 所示?？梢园l(fā)現(xiàn)，曲線的變化趨勢與第1 種方案（彈性模量分?jǐn)?shù)比理論）得到的結(jié)果類似，即AL先增大后減小，5 GPa 處出現(xiàn)轉(zhuǎn)折。整體上看，高壓下HBT 晶體彈性模量的各向異性程度減弱。由此可知，AC、AU和AL三者在描述HBT 晶體各向異性程度上是等價(jià)的。

圖 8 AU 和AC 的各向異性關(guān)系Fig. 8 Anisotropy diagram of the AU and AC

圖 9 AL 隨壓強(qiáng)的變化關(guān)系Fig. 9 Pressure as a function of anisotropy index AL

從圖7、圖9 可以明顯看出，5、25 和55 GPa 下曲線有較為明顯的變化。為了進(jìn)一步考察這3 個(gè)壓強(qiáng)下HBT 晶體模量的各向異性程度，利用ELATE 開源軟件[31]，將這3 個(gè)壓強(qiáng)下的彈性系數(shù)6×6 矩陣元36 個(gè)數(shù)值全部代入，繪制出剪切模量的二維分布，如圖10 所示。因單位換算，圖10 中坐標(biāo)軸均被放大了10 倍。圖10(a)、圖10(b)、圖10(c)分別顯示了5、25 和55 GPa 下剪切模量的二維分布。每個(gè)壓強(qiáng)下的3 幅圖分別代表剪切模量三維分布在xy、xz、yz3 個(gè)面上的投影。藍(lán)線和綠線分別代表兩個(gè)方向剪切模量的矢量投影。從圖10 可以看出：壓強(qiáng)越大，每個(gè)面上投影矢量的各向異性越強(qiáng)；在確定的壓強(qiáng)下，每個(gè)面上的投影也表現(xiàn)出差異性。綜上所述，HBT 晶體剪切模量的各向異性受壓強(qiáng)的影響顯著。

圖 10 不同壓強(qiáng)下剪切模量各向異性的二維投影Fig. 10 Anisotropic two-dimensional projection of shear modulus at different pressures

3 結(jié) 論

采用基于密度泛函理論的第一性原理方法計(jì)算了HBT 晶體在零溫零壓及高壓下的結(jié)構(gòu)，研究了高壓下HBT 的力學(xué)穩(wěn)定性和彈性性質(zhì)，具體分析了HBT 晶體的彈性常數(shù)Cij、彈性模量（B、G、E）、聲速（ vp、 vs、 vm）和德拜溫度Θ 隨壓強(qiáng)的變化關(guān)系。在此基礎(chǔ)上，基于3 種不同理論研究了彈性模量的各向異性。結(jié)果表明，彈性常數(shù)和模量都隨壓強(qiáng)的增加而增加。根據(jù)B/G 估測了HBT 晶體在高壓下的韌脆性。彈性模量各向異性參量隨壓強(qiáng)表現(xiàn)出不同的變化規(guī)律，總體而言，高壓下HBT 晶體的彈性模量各向異性程度減弱。