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

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

隨著化石燃料的減少與溫室效應(yīng)的加劇，能源問題日益嚴(yán)重。使用相變材料（PCM）來存儲(chǔ)清潔能源是當(dāng)今社會(huì)的熱門話題。金屬和合金由于其導(dǎo)熱系數(shù)高、儲(chǔ)熱密度高、性能穩(wěn)定等優(yōu)點(diǎn)得到了學(xué)者們的青睞，為解決儲(chǔ)能問題提供了的重要思路[1-2]。

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

1 模型與方法

1.1 模型構(gòu)建與優(yōu)化

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

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

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

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

1.2 熱物性計(jì)算

分子動(dòng)力學(xué)計(jì)算密度公式，其表達(dá)見式（3）～式（4）

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

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

比熱容的計(jì)算是通過系統(tǒng)總能量對溫度的求導(dǎo)并除以總質(zhì)量得到的，公式表達(dá)如式（5）所示[13]

式中，m是系統(tǒng)總質(zhì)量，是總能量對溫度的求導(dǎo)。

計(jì)算熱導(dǎo)率的方法包括平衡分子動(dòng)力學(xué)模擬（EMD）和非平衡分子動(dòng)力學(xué)模擬（NEMD）。在本文中采用非平衡分子動(dòng)力學(xué)模擬計(jì)算熱導(dǎo)率。首先建立溫度梯度并等待足夠長的時(shí)間以使系統(tǒng)達(dá)到穩(wěn)定狀態(tài)。達(dá)到穩(wěn)定狀態(tài)后，測量系統(tǒng)的溫度梯度并計(jì)算熱導(dǎo)率。恒溫力由式（6）給出[14]

式中，mi表示質(zhì)量，k(ri)表示所處的相對位置，F(xiàn)Гk(ri)表示進(jìn)入體系的能量，KГk(ri)和VГk(ri)表示非轉(zhuǎn)變動(dòng)能和質(zhì)心速度。

產(chǎn)生溫度梯度的方法如下：采用周期性邊界條件，在X方向上被均分為20個(gè)塊，第1塊作為熱源，第11塊作為熱匯。在系統(tǒng)達(dá)到某一溫度平衡之后，移除整個(gè)熱浴，分別對熱源和熱匯區(qū)域施加局部Nose-Hoover熱浴，并施加適當(dāng)功率。在系統(tǒng)達(dá)到穩(wěn)定狀態(tài)后，測量溫度梯度，得到熱導(dǎo)率。

雙體分布函數(shù)（pair correlation function）是分子動(dòng)力學(xué)中研究液體、晶體和非晶體的重要分析方法，是理論與實(shí)際聯(lián)系的重要工具。其定義為

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

2 結(jié)果與討論

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

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

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

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

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

從圖5可以看出，Mg-Cu合金的比熱容在熔化時(shí)突然變大，這表示此時(shí)發(fā)生了熔化，從而確定了其熔點(diǎn)。在熔化完成后，比熱容又突然減小到熔化前的大小。在相變前后，比熱容隨溫度變化不大，保持平穩(wěn)狀態(tài)。隨著Cu含量增加，Mg-Cu合金的熔點(diǎn)先減小后增加，Mg-31Cu合金的熔點(diǎn)最低。隨著Cu含量增加，其比熱容呈降低的趨勢。這是因?yàn)樵诔叵陆饘冁V的比熱容為1.02J?g-1?K-1，金屬銅的比熱容為0.39 J?g-1?K-1，鎂的比熱容要大于銅的，而合金的比熱容可以作為各組分比熱容與其質(zhì)量分?jǐn)?shù)的乘積之和（Neumann-Kopp定律），所以 Mg含量越高，合金的比熱容越大。

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

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

圖6 熱導(dǎo)率隨溫度變化Fig.6 Thermal conductivity variations with temperature

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

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

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

3 結(jié) 論

本文采用EAM勢分子動(dòng)力學(xué)方法模擬了Mg-Cu合金的熱物性質(zhì)，分析了其微觀結(jié)構(gòu)。主要結(jié)論如下。

圖7 Mg-Cu合金的雙體分布函數(shù)（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含量的增加，其密度單調(diào)增加，與實(shí)驗(yàn)值一致。根據(jù)能量-溫度曲線分析、比熱容-溫度曲線得到，隨著Cu含量的增加，熔點(diǎn)先降低后升高，Mg-31Cu合金的熔點(diǎn)最低。在相變前后，比熱容隨溫度變化不大，但比熱容隨著Cu含量增加呈降低的趨勢。其熔化焓隨著Cu含量的增加，先增加后減小。Mg-31Cu合金的熔化焓最大，為156.46J/g。熱導(dǎo)率隨Cu含量的增加而降低，但即使Cu含量達(dá)到45%，其熱導(dǎo)率仍高達(dá)100.16 W?(m?K)-1。

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