李 佐，張鳳玲，廖大麟

（貴州工程應(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%的氮元素，爆炸燃燒后釋放的氣體主要為氮氣、二氧化碳和水，因此被稱為“綠色”含能材料[8-11]。本研究采用基于第一性原理的準諧振近似方法，計算固態(tài)HBT 晶體在0～120 GPa 高壓下的彈性性質(zhì)，同時分析其彈性的各向異性特征。

1 理論與計算方法

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

采用基于第一性原理的Quantum ESPRESSO 軟件包中的PWSCF 模塊[12]計算電子結(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]。為增加計算精度，平面波基函數(shù)的截斷能取為544 eV，布里淵空間采用Monkhorst-Pack 方法[17]，k 點方案為6 × 6 × 3，總能自洽，離子位移最大值和應(yīng)力最大值都能達到設(shè)置的收斂精度。在高壓下幾何結(jié)構(gòu)優(yōu)化的基礎(chǔ)上，利用準諧振近似方法計算HBT 晶體的彈性性質(zhì)，詳細的計算過程在軟件Thermo_pw 中體現(xiàn)[18]。

1.2 彈性常數(shù)

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

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

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

表 1 HBT 晶體晶格參數(shù)的計算值和實驗值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 理論和實驗得到的壓強-體積比關(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ǔ)上計算彈性系數(shù)和彈性模量。三斜結(jié)構(gòu)的彈性常數(shù)矩陣包含21 個獨立矩陣元。Born[22]d 給出了晶體材料力學(xué)穩(wěn)定性的充分必要條件為彈性常數(shù)矩陣C 正定。按照線性代數(shù)知識可知，對稱矩陣C 正定的充要條件為C 的各階順序主子式都為正，即

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

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

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

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

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

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

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

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

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

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

圖 6 不同壓強下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]提出的體彈模量、楊氏模量和剪切模量分數(shù)比各向異性的概念，定義為

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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