Mg-Cu合金熱物性的分子動力學計算

余嘉鵬，程曉敏,2，李元元，李 蓓,徐 虹

（1武漢理工大學材料科學與工程學院，湖北 武漢 430070；2黃岡師范學院機電與汽車工程學院，湖北 黃岡 438000；3洛林大學LCP-A2MC實驗室,法國 梅斯 57070）

隨著化石燃料的減少與溫室效應的加劇，能源問題日益嚴重。使用相變材料（PCM）來存儲清潔能源是當今社會的熱門話題。金屬和合金由于其導熱系數高、儲熱密度高、性能穩定等優點得到了學者們的青睞，為解決儲能問題提供了的重要思路[1-2]。

目前，針對能源領域已經研究了許多合金系統，常見的有Mg，Al，Cu，Si，Zn等[3]。其中，鎂基合金由于其優異的熱物性能，國內外許多學者對其進行了許多的研究。P BLANCO-RODRíGUEZ和J RODRíGUEZ-ASEGUINOLAZA等[4-6]研究了Mg基二元三元合金儲熱材料的熱性能。程曉敏等[7]研究了Mg-Cu-Zn相變儲熱材料，并用EET理論對其儲熱機理提出了解釋。孫正等[8]選取了7種元素加入Mg基體中，得到了多種二元和三元高溫相變儲熱合金。在上述文獻中， Mg基合金儲熱材料通過實驗方法對其熱物性能進行了廣泛研究。本文通過分子動力學的方法構建了不同組分的Mg-Cu合金塊體，對其熔化焓、熱導率等熱物性參數進行了模擬，并且用分子動力學方法表征了Mg-Cu合金的微觀結構，分析Mg-Cu合金熱物性能變化的微觀機理。為儲熱材料的研究提供了一種新思路。

1 模型與方法

1.1 模型構建與優化

晶體鎂具有密排六方結構，通過晶體鎂可以構建出不同粒徑的鎂塊體，如圖1(a)所示。隨后按照質量比重，將晶體鎂中的原子隨機替換銅原子，如圖1(b)所示。本文模擬構建了不同組分Mg-Cu合金塊體，采用共軛梯度算法對得到的Mg-Cu合金塊體團簇進行能量最小化，使團簇能量達到最小狀態。所有模擬均使用LAMMPS[9]軟件進行，OVITO[10]用于后期處理。

本模擬采用嵌入原子勢（EAM）[11]，它能很好地描述金屬及其合金之間的相互作用，得到金屬的力學、熱力學性能。本文采用了Howard Sheng獲得的Mg-Cu合金的勢函數，該系統的總能量具有以下公式[12]：

式中，ρh,i是由原子i與系統中剩余原子所產生的主電子密度；Fi(ρh,i)代表原子i嵌入時產生的嵌入能；?ijRij代表原子i與原子j之間的作用勢；Rij是它們之間的距離。

式中，aj(R)是由原子j貢獻的電子密度。對整個合金體系，采用周期性邊界條件，在NPT系綜下進行模擬。進行牛頓運動方程求解時采用Verlet積分算法，時間步長設置為1fs，溫度阻尼系數設定為時間步長的100倍。設置目標溫度為300K，并弛豫足夠長的時間以達到熱力學平衡，使結構更穩定。圖2所示為弛豫過程中溫度和能量的變化，可以看到500ps后曲線趨于平衡，整個體系達到了熱力學平衡,得到了穩定的合金模型。

1.2 熱物性計算

分子動力學計算密度公式，其表達見式（3）～式（4）

式中m為總質量，N1為模型中Mg原子數，N2為模型中Cu原子數，M1為Mg的摩爾質量，M2為Cu的摩爾質量，V為計算給出的總體積。為了避免統計誤差，對于每次得到的密度取時間平均。

通過能量-溫度曲線與比熱容-溫度曲線兩條曲線來計算Mg-Cu合金的熔點。相變發生點附近，系統的總能量會發生突變，這就說明體系發生了相變。但是通過系統總能量計算出的熔點不夠準確，通過計算比熱容能得到更精確的熔點。熔點的模擬是采用NPT系綜，升溫速率為1×1012K/s。

比熱容的計算是通過系統總能量對溫度的求導并除以總質量得到的，公式表達如式（5）所示[13]

式中，m是系統總質量，是總能量對溫度的求導。

計算熱導率的方法包括平衡分子動力學模擬（EMD）和非平衡分子動力學模擬（NEMD）。在本文中采用非平衡分子動力學模擬計算熱導率。首先建立溫度梯度并等待足夠長的時間以使系統達到穩定狀態。達到穩定狀態后，測量系統的溫度梯度并計算熱導率。恒溫力由式（6）給出[14]

式中，mi表示質量，k(ri)表示所處的相對位置，FГk(ri)表示進入體系的能量，KГk(ri)和VГk(ri)表示非轉變動能和質心速度。

產生溫度梯度的方法如下：采用周期性邊界條件，在X方向上被均分為20個塊，第1塊作為熱源，第11塊作為熱匯。在系統達到某一溫度平衡之后，移除整個熱浴，分別對熱源和熱匯區域施加局部Nose-Hoover熱浴，并施加適當功率。在系統達到穩定狀態后，測量溫度梯度，得到熱導率。

雙體分布函數（pair correlation function）是分子動力學中研究液體、晶體和非晶體的重要分析方法，是理論與實際聯系的重要工具。其定義為

式中，i與j表示不同的原子，Ni與Nj分別表示兩種原子數，gi,j(r)表示以α原子為中心，距離在r→r+Δr內找到另一個原子的概率，L表示元胞的邊長，Δr表示計算步長，nαj(r)表示以i原子為中心，距離在r→r+Δr內發現j原子的數量。

2 結果與討論

圖3所示為Mg-Cu合金的密度隨溫度的變化。初始溫度為300 K，之后逐漸升溫到700 K。在此過程中，Mg-Cu合金的密度單調遞減，變化趨勢平緩。隨著Cu含量的增加，密度單調遞增。

圖3 密度隨溫度變化Fig.3 Density variations with temperature

下文對Mg-Cu合金的相變行為進行研究，模擬了不同組分下，其相變溫度區間、熔點的變化情況。如圖4所示為Mg-Cu合金熔化過程中總能量隨溫度的變化。在升溫熔化的過程中，系統總能量在熔點處發生了明顯的跳躍，即觀察到了Mg-Cu合金相變的過程，在相變前后，系統總能量均呈線性增長。從圖4可以得到Mg-Cu合金發生相變的溫度范圍，隨著Cu含量的增加，Mg-Cu合金的熔點先增大后減小。

為了得到Mg-Cu合金更準確的熔點，本文模擬了相變過程中Mg-Cu合金比熱容的變化情況。圖5表示了不同組分的Mg-Cu合金塊體比熱容隨溫度的變化。

圖4 熔化過程中總能量的溫度曲線Fig.4 Total energy changing in process of melting

從圖5可以看出，Mg-Cu合金的比熱容在熔化時突然變大，這表示此時發生了熔化，從而確定了其熔點。在熔化完成后，比熱容又突然減小到熔化前的大小。在相變前后，比熱容隨溫度變化不大，保持平穩狀態。隨著Cu含量增加，Mg-Cu合金的熔點先減小后增加，Mg-31Cu合金的熔點最低。隨著Cu含量增加，其比熱容呈降低的趨勢。這是因為在常溫下金屬鎂的比熱容為1.02J?g-1?K-1，金屬銅的比熱容為0.39 J?g-1?K-1，鎂的比熱容要大于銅的，而合金的比熱容可以作為各組分比熱容與其質量分數的乘積之和（Neumann-Kopp定律），所以 Mg含量越高，合金的比熱容越大。

下文對Mg-Cu合金的熱導率進行了模擬研究。如圖6所示，隨著溫度的上升，Mg-Cu合金的熱導率呈單調增加的趨勢。隨著Cu含量的增加，Mg-Cu合金的熱導率單調降低。熱導率降低有兩個原因：一方面，隨著Cu含量的增加，第二相中Mg2Cu的含量增加，形成了許多新的界面和晶格的畸變，使其成為電子散射的來源，導致合金的熱導率降低；另一方面；Cu原子的加入破壞了Mg原子的周期性排列，由于Cu和Mg原子的半徑不同，引起了晶格畸變，增強電子的散射，阻礙了自由電子在晶格中的運動，從而導致了合金的熱導率降低[15]。

圖5 比熱容隨溫度變化Fig.5 CP variations with temperature

圖6 熱導率隨溫度變化Fig.6 Thermal conductivity variations with temperature

表1 Mg-Cu合金的熱物性 Table1 Thermal properties of Mg-Cu alloys

表1給出了模擬得到Mg-Cu合金的熔點、熔化焓、比熱容和熱導率的具體值?？梢钥吹?，在這7組合金中，隨著Cu含量的增加，Mg-Cu合金的熔化焓先增大后減小，Mg-31Cu合金的熔化焓最大，為156.46 J/g。而金屬鎂的熔化焓是要大于金屬銅的，如果按照質量比重計算熔化焓，則隨著Cu含量的增加，其熔化焓變化應該呈現單調遞減的趨勢。我們認為熔化焓出現這種現象的原因可能是Mg-Cu合金在相變中生成了α+Mg2Cu共晶組織，α+Mg2Cu共晶組織的含量可能是影響了Mg-Cu合金系統的熔化焓的一個主要原因[16]。

下文用雙體分布函數的表征方法對合金的微觀結構進行模擬分析。圖7（a）、（b）、（c）分別表示了在室溫下Mg-24Cu、Mg-31Cu和Mg-40Cu三種合金的雙體分布函數曲線。我們對比三種合金的gMg-Cu(r)曲線發現，Mg-31Cu與Mg-40Cu合金的第一峰比Mg-24Cu合金的要更高，峰位也向右偏移，這說明Mg-31Cu與Mg-40Cu合金中Mg原子與Cu原子成鍵數比Mg-24Cu合金的要多。對比三種合金的gMg-Mg(r)曲線，Mg-24Cu合金的第一峰比Mg-31Cu與Mg-40Cu合金要高。對比三種合金的gCu-Cu(r)曲線，Mg-40Cu合金的第一峰遠比Mg-24Cu與Mg-31Cu合金要高，這說明Mg-40Cu合金中Cu-Cu成鍵最多。結合Mg-Cu二元合金相圖，Mg-24Cu合金與Mg-31Cu合金生成的相為初生的α-Mg相與α-Mg+Mg2Cu共晶組織，Mg-31Cu合金中的α-Mg+Mg2Cu共晶組織要多于Mg-24Cu合金。而Mg-40Cu合金生成的相為α-Mg+Mg2Cu共晶組織和Mg2Cu相。

3 結 論

本文采用EAM勢分子動力學方法模擬了Mg-Cu合金的熱物性質，分析了其微觀結構。主要結論如下。

圖7 Mg-Cu合金的雙體分布函數（a）Mg-24Cu;（b）Mg-31Cu;（c）Mg-40CuFig.7 Pair distribution function of Mg-Cu alloys：(a) Mg-24Cu;(b) Mg-31Cu;(c) Mg-40Cu

（1）Mg-Cu合金的密度隨溫度的升高而降低，隨著Cu含量的增加，其密度單調增加，與實驗值一致。根據能量-溫度曲線分析、比熱容-溫度曲線得到，隨著Cu含量的增加，熔點先降低后升高，Mg-31Cu合金的熔點最低。在相變前后，比熱容隨溫度變化不大，但比熱容隨著Cu含量增加呈降低的趨勢。其熔化焓隨著Cu含量的增加，先增加后減小。Mg-31Cu合金的熔化焓最大，為156.46J/g。熱導率隨Cu含量的增加而降低，但即使Cu含量達到45%，其熱導率仍高達100.16 W?(m?K)-1。

（2）表征了Mg-Cu合金在常溫下的微觀結構，結果表明，Mg-Cu合金中不同原子之間存在更強的相互作用；結合Mg-Cu二元合金相圖分析得到，Mg-24Cu合金與Mg-31Cu合金生成的相為初生的α-Mg相與α-Mg+Mg2Cu共晶組織，Mg-31Cu合金中的α-Mg+Mg2Cu共晶組織要多于Mg-24Cu合金；Mg-40Cu合金生成的相為α-Mg+Mg2Cu共晶組織和Mg2Cu相；α-Mg+Mg2Cu共晶組織是導致Mg-31Cu合金的熔化焓的增加以及熱導率變化的主要原因。