張 玥， 劉 泉

(安徽大學 物理與光電工程學院，安徽 合肥 230601)

體積彈性模量與二階彈性常數(shù)是表征熱彈性性質(zhì)的重要物理量之一，要想全面了解地球內(nèi)部的成分和熱狀態(tài)，就需要了解所有備選礦物在極端壓強和溫度條件下的彈性特性.長期以來，人們建立了許多計算固體彈性常數(shù)的理論模型，并對其與溫度和壓強的關(guān)系進行了深入的研究[1-4].

2010年，基于彈性常數(shù)對溫度的假設(shè)性依賴關(guān)系，Wang等人提出了一種計算有限溫度下彈性剛度系數(shù)的第一性原理方法，結(jié)果表明預測值與已有的實驗測量值一致[5].同年，Zen等人使用第一性原理方法研究了Ni，Ti和NiTi合金的結(jié)構(gòu)和彈性性能，其中彈性常數(shù)是通過使用密度泛函技術(shù)評估能量與小應變的關(guān)系來計算的[6].2013年，Ono等人采用從頭算分子動力學方法研究了立方體鈣鈦礦在高壓和高溫下的彈性，計算結(jié)果表明立方體鈣鈦礦是地球下地幔的一種次要礦物[7]. 2018年，Liu等人利用基于密度泛函微擾理論的第一性原理方法，研究了直至2500 K、100 GPa高溫、高壓下黃鐵礦的等溫體積模量、熱膨脹系數(shù)、熱容和熵等熱力學性質(zhì),與實驗結(jié)果保持了很好的一致性[8].2019年，Ulian以及Valdrè采用從頭算量子力學方法計算了2種礦物的狀態(tài)方程(EoS)和二階彈性常數(shù)，理論結(jié)果與現(xiàn)有文獻中實驗觀察到的一般趨勢一致，并進一步擴展了對兩相力學性能的認識[9]. 2020年，Argaman與Makov用第一性原理各向異性準調(diào)和方法計算了HCP鈦的彈性常數(shù)，其中包括聲子譜作為應變的函數(shù)計算，除C13外，彈性常數(shù)的溫度依賴性與實驗結(jié)果一致[10]. 2021年，Yang等人采用第一性原理計算，研究Mg，Be，Ti，Zn，Zr和Cd的三階彈性常數(shù)和機械性能，結(jié)果與之前的計算和實驗值更加吻合[11].

在這里，我們應該指出的是，無論是第一性原理計算還是從頭算分子動力學方法[5-11]，其間均包含大量繁雜的計算工作.尋找一種簡單直接的方法計算多種礦物質(zhì)在不同溫度下的體積彈性模量及二階彈性常數(shù)，成為當今地球物理學中迫切需要解決的問題.

本文中，我們通過細致分析大量地幔礦物質(zhì)高溫下的熱力學參量，基于泰勒級數(shù)展開獲得熱膨脹系數(shù)與溫度之間的唯象依賴性關(guān)系式，并結(jié)合熱力學公式與Tallon模型[12]，計算了MgAl2O4，Mg2SiO4及Al2O3這3種典型地幔礦物質(zhì)在不同溫度下的體積彈性模量及二階彈性常數(shù).

1 分析方法

為了探尋體積彈性模量(BT)對溫度T的依賴性，Anderson引入了著名的Anderson-Grüneisen參數(shù)δT，其定義如下[1]：

(1)

其中：P為壓強，α為體積熱膨脹系數(shù).方便起見，我們將δT寫為δ，將BT寫成B.根據(jù)α的定義(其中V表示體積)：

(1)式可以寫為

(2)

Chopelas 和 Boehler[13]提出了一個等溫Anderson-Grüneisen參數(shù)δ與體積V之間的線性依賴關(guān)系式：

(3)

對給定的晶體，A是一個常數(shù).結(jié)合(2)～(3)式并考慮熱壓項，Kumar得到下列表達式[14-15]：

[1-(δ0+1)α0(T-T0)]，

(4)

式中α0、B0及δ0分別為α、B和δ在初始溫度T=T0=300 K時的值.Tallon 已經(jīng)發(fā)現(xiàn)(4)式對于任何一個彈性模量都成立[9]，從而可以用M替代B得到

式中M代表任何一個彈性常數(shù)，如C11,C12,C44或者B.從以上分析可知，要求得任意溫度下的彈性常數(shù)，必須首先求得體積彈性模量與溫度之間的關(guān)系式.以下分析將幫助我們尋找一種合適的B-T關(guān)系式.

通過分析不同溫度下大量礦物質(zhì)的熱膨脹系數(shù)實驗數(shù)據(jù)，基于泰勒級數(shù)展開式，將熱膨脹系數(shù)與溫度之間的依賴性關(guān)系表示為

α=α0+α′0(T-T0)+

(5)

所以(5)式可以改寫為

(6)

將(6)式代入(1)式并積分，積分時當溫度從T0到T時，相應的體積彈性模量從B0到B.最后得到B與T之間的關(guān)系式：

(7)

根據(jù)Tallon的推廣方法[12]可以將(7)式推廣為

(8)

(9)

表1 計算所需輸入數(shù)據(jù)(取自文獻[1])

2 計算結(jié)果和討論

(5)式已被Kumar[14-15]推廣用于計算高溫下大量地幔礦物質(zhì)的彈性常數(shù).在本文中，我們首先運用(5)式和(7)式計算不同溫度下NaCl晶體的體積彈性模量.計算結(jié)果和實驗數(shù)據(jù)見表2,實驗數(shù)據(jù)取自文獻[2]，其中300 K時的Anderson-Grüneisen參數(shù)根據(jù)文獻[2]取作5.95.由表2可知，與實驗數(shù)據(jù)相比，利用(7)式計算的結(jié)果誤差更小，即(7)式更準確地反映了體積彈性模量隨溫度的變化規(guī)律.

表 2 分別利用 (5)式和 (7) 式計算得到的NaCl體積彈性模量隨溫度的變化值

對于MgAl2O4，存在3個二階彈性常數(shù)：C11,C12和C44.從圖1可以看出，彈性常數(shù)計算值與實驗值吻合度較高，最大誤差僅為0.3%.此外，相對于C12和C44，C11隨溫度變化的幅度更大.這是因為C11代表縱向彈性，而縱向的應變只引起體積改變而無形狀改變，體積的改變又與溫度變化密切相關(guān)，從而導致了C11的巨大變化.另一方面，C12和C44代表橫向彈性即剪切常數(shù)，橫向應變只引起形狀變化而無體積變化，所以C12和C44對溫度并不敏感.

圖1 不同溫度下MgAl2O4的體積彈性模量(GPa)及二階彈性常數(shù)

(9)式還可用以計算更復雜的礦物質(zhì)，如Mg2SiO4及Al2O3.對于Mg2SiO4，存在9個二階彈性常數(shù)，而Al2O3存在6個二階彈性常數(shù).從圖2與圖3可以看出，根據(jù)(7)式與(9)式計算得到的體積彈性模量和二階彈性常數(shù)與實驗數(shù)值的一般趨勢也保持一致.當溫度不超過1 000 K時，計算值與實驗值完全一致.隨著溫度的繼續(xù)升高，少數(shù)二階彈性常數(shù)的誤差有所增大.Mg2SiO4最大的誤差出現(xiàn)在1 700 K時的C33，為3.9%；而Al2O3的最大誤差出現(xiàn)在1 800 K時的C33，為3.8%.

圖2 不同溫度下Mg2SiO4的體積彈性模量(GPa)及二階彈性常數(shù)

3 結(jié)語

本文通過分析大量地幔礦物質(zhì)高溫下的實驗數(shù)據(jù)，結(jié)合熱力學關(guān)系式及高溫下礦物質(zhì)熱膨脹系數(shù)的唯象公式，提出了一種簡單而直接的數(shù)值方法以計算不同溫度下固體彈性常量.同時，本文所用方法不依賴于礦物質(zhì)的結(jié)構(gòu)，因此還可應用于更復雜的地幔礦物質(zhì).