NaNO3-KNO3-NaNO2三元混合相變?nèi)埯}結(jié)構(gòu)與物性的分子動力學(xué)模擬

倪海歐，孫 澤，路貴民，于建國

?

NaNO3-KNO3-NaNO2三元混合相變?nèi)埯}結(jié)構(gòu)與物性的分子動力學(xué)模擬

倪海歐，孫 澤，路貴民，于建國

（華東理工大學(xué)國家鹽湖資源綜合利用工程技術(shù)研究中心，上海200237）

相變?nèi)埯}目前是太陽能熱發(fā)電的重要傳蓄熱介質(zhì)，其結(jié)構(gòu)和性質(zhì)之間的定量關(guān)系是相變?nèi)埯}研究的熱點。本文采用分子動力學(xué)方法，利用Buckingham勢函數(shù)，對NaNO3-KNO3-NaNO2三元混合熔鹽的結(jié)構(gòu)和性質(zhì)進行了研究。計算了混合熔鹽在不同溫度下的徑向分布函數(shù)、配位數(shù)、角分布函數(shù)等結(jié)構(gòu)信息，以及密度、剪切黏度、熱導(dǎo)率、比熱容等性質(zhì)。結(jié)果表明，計算得到的各項性質(zhì)與文獻值吻合較好，驗證了勢函數(shù)和計算方法的可靠性，為計算機模擬指導(dǎo)熔鹽配方、熔鹽構(gòu)效關(guān)系定量研究奠定了基礎(chǔ)。

分子動力學(xué)；熔鹽結(jié)構(gòu)；熔鹽物性

隨著能源問題越來越成為當(dāng)今世界所關(guān)注的一大重要議題，太陽能光熱發(fā)電作為一種綠色可持續(xù)的能源形式，正受到越來越廣泛的關(guān)注。我國已在“十二五”期間啟動示范項目，“十三五”規(guī)劃到2020年光熱發(fā)電裝機容量達到500萬千瓦[1]，光熱發(fā)電正進入一個快速發(fā)展期。

熔鹽作為一種廉價、高效的傳蓄熱介質(zhì)，在太陽能熱發(fā)電等領(lǐng)域有著重要的應(yīng)用[2-3]。但是，由于熔鹽的高溫特性，其物性的實驗測試成本高、難度大、精度低，對測試儀器的傷害也很大。而隨著計算機科學(xué)的飛速發(fā)展，計算機模擬可以處理的體系越來越大，計算精度也越來越高，在很多情況下可以作為傳統(tǒng)實驗的替代和補充，尤其適合熔鹽這種高溫高腐蝕的場合。LYNDEN-BELL等[4]用分子動力學(xué)方法計算了NaNO2的晶格振動。PENG等[5]用第一性原理方法計算了KNO2固體的Debye溫度、比熱容、聲子速度、聲子平均自由程、聲子熱導(dǎo)率等性質(zhì)。JAYARAMAN等[6]利用分子動力學(xué)方法計 算了Li、Na、K三種堿金屬硝酸鹽的密度、黏度、 熱導(dǎo)率、比熱容、熔點等性質(zhì)。URAHATA等[7]采用帶有極化模型的分子動力學(xué)方法計算了LiNO3的一些性質(zhì)，并比較了極化模型與經(jīng)典模型計算結(jié)果的差異。

本文通過經(jīng)典分子動力學(xué)模擬的方法計算了質(zhì)量比為7∶53∶40的NaNO3-KNO3-NaNO2三元混合熔鹽。這種熔鹽通常被稱為HITEC，在工業(yè)上有著廣泛的應(yīng)用。本文計算了HITEC熔鹽的結(jié)構(gòu)與性質(zhì)，包括徑向分布函數(shù)、角分布函數(shù)、密度、剪切黏度、熱導(dǎo)率、比熱容等，并與文獻中的實驗值進行了比較，驗證了勢函數(shù)和計算方法的可靠性。

1 分子動力學(xué)計算

1.1 力場模型

采用Buckingham勢函數(shù)來描述模擬體系分子間庫侖相互作用、短程排斥作用和范德華耗散作用，如式(1)所示

式中，q、q為原子電荷，為原子間距離，r為截斷半徑，、為短程斥力的力場參數(shù)，為耗散系數(shù)?；旌先埯}中，硝酸鈉和硝酸鉀采用文獻[8]中的勢函數(shù)，亞硝酸鈉采用文獻[4]中的勢函數(shù)。勢函數(shù)參數(shù)見表1（硝酸根中的N和O以N3、O3表示，亞硝酸根中的以N2、O2表示，下文均用這種方法表示）。對于交叉項，采用如下交叉法則：A= (AA)1/2，C= (CC)1/2，2/ρ=1/ρ+1/ρ。

硝酸根與亞硝酸根的分子內(nèi)作用力的力常數(shù)由拉曼光譜數(shù)據(jù)計算得出[9-10]，N—O鍵采用Harmonic簡諧勢，如式(2)所示，式中為N—O鍵長，0為平衡鍵長，K為力常數(shù)；O—N—O鍵角采用Charmm勢，如式(3)所示，其中第一項為O—N—O鍵角彎曲的簡諧項，第二項為O—O非鍵排斥項。式中為O—N—O鍵角，0為平衡鍵角，K為力常數(shù)，為O—O距離，UB為平衡距離，UB為Urey Bradley力常數(shù)。具體參數(shù)見表2。

（3）

1.2 計算方法

模擬計算使用開源計算軟件Lammps進行。采用周期性邊界條件，總粒子數(shù)為5540個，其中含有685個Na+、543個K+、628個、600個。勢函數(shù)截斷半徑取模擬盒子邊長的一半。長程作用力采用pppm求和法，計算精度設(shè)置為10-4。時間步長取1fs。

表1 分子間勢函數(shù)參數(shù)

表2 分子內(nèi)勢函數(shù)參數(shù)

2 結(jié)果與討論

2.1 徑向分布函數(shù)

徑向分布函數(shù)（radical distribution function，RDF）是描述流體結(jié)構(gòu)特征最重要的函數(shù)，其定義如式(4)所示，其中β為β粒子的數(shù)密度，αβ()表示位于以α粒子為中心、為半徑的球體內(nèi)β粒子的平均數(shù)。

圖1是陽-陽離子對的RDF。從圖中可以看出，隨著溫度的上升，RDF第一峰的高度逐漸降低，第一峰峰谷的高度逐漸增大，峰變得更平，峰位置右移，而RDF的整體形貌沒有發(fā)生明顯變化。圖2是陽-陰離子對的RDF，同樣的，隨著溫度的升高，峰趨向于展平，但峰位置沒有明顯變化。這是由于溫度的升高，離子之間相互作用減弱，分子運動隨機性變大，分子微觀團簇變得有所松散。而溫度的升高使得同號離子之間距離變大，導(dǎo)致峰位置右移，但異號離子之間庫侖引力較強，溫度的升高對其離子間距離并沒有明顯影響。

圖3是623 K下異號離子對的RDF。從圖中可以看出，Na-N3和K-N3的峰高基本一致，而Na-N3峰位置比K-N3更小，兩者峰位置之差是由于Na的離子半徑小于K，導(dǎo)致Na-N之間距離更近。而Na-N2的峰高明顯高于K-N2，兩者第一峰位置的差異也相對較大，這可能是源于Na+與的相互作用更強，兩者結(jié)合得更加緊密。

配位數(shù)是描述熔鹽結(jié)構(gòu)最重要得參數(shù)之一，如式(5)所示，可由徑向分布函數(shù)得積分得到。表3列出了M-N離子對的配位數(shù)。可以看出，隨著溫度的升高，所有的離子對的配位數(shù)均變小了，這說明隨著溫度的升高，微觀層面上離子團簇變得松散。同時，M-N3離子對的配位數(shù)比M-N2離子對大，這是由于體系中的數(shù)量多于。對不同的陽離子，Na-N的配位數(shù)小于K-N，而在徑向分布函數(shù)中Na-N的峰高于K-N，這說明雖然K的第一配位殼層內(nèi)陰離子個數(shù)比Na多，但是和Na相比結(jié)合的更為松散。

表3 陽-陰離子的配位數(shù)

2.2 角分布函數(shù)

徑向分布函數(shù)在描述液體局域結(jié)構(gòu)中起著至關(guān)重要的作用，但是其只描述了粒子之間的配對狀況，而沒有描述其團簇結(jié)構(gòu)內(nèi)部粒子之間的取向信息。為了獲取團簇內(nèi)部的取向信息，在徑向分布函數(shù)的基礎(chǔ)上，本文還計算了體系陰-陽-陰離子的角度分布函數(shù)（angular distribution function，ADF），以反映“鍵”的取向[11]。

圖4為混合熔鹽中N-M-N“鍵角”的角分布函數(shù)（M指代Na和K）。從圖中可以看出N-M-N“鍵角”主要分布在40°～180°之間，在90°附近分布最為密集，說明在熔鹽中M-N“鍵角”的取向并不是無規(guī)律的，而是傾向于形成松散的八面體配位結(jié)構(gòu)。

圖5選取了623 K下N-M-N角分布函數(shù)，可以看出對相同的陽離子，N3-M-N3和N2-M-N2的角分布沒有明顯差別。而對不同的陽離子，N-Na-N的峰比N-K-N的峰更高，且更接近90°，說明Na-N的八面體配位構(gòu)型比K-N更為顯著。N-K-N的峰小于90°也和K-N的配位數(shù)超過6吻合。

2.3 密 度

密度計算在NPT系綜下進行，每個溫度點共計算200萬步，然后對每一步輸出的密度值進行平均，得到該溫度下熔鹽密度的模擬值。模擬值與文獻值的對比如圖6所示，兩者的偏差小于3%。密度的文獻值采用Serrano-López[12]總結(jié)的關(guān)聯(lián)式，以下剪切黏度、比熱容的文獻值均采用該文獻的值。

2.4 剪切黏度

熔鹽的剪切黏度采用平衡分子動力學(xué)方法計算，根據(jù)Green-Kubo公式，由模擬體系應(yīng)力張量自相關(guān)函數(shù)（Stress Tensor Auto-Correlation Function）計算出體系的剪切黏度，如式(6)所示。式中B為波爾茲曼常數(shù)，、為體系的溫度和體積，S為應(yīng)力張量的分量。應(yīng)力張量的每一個非對角分量（S，S，S）都可以計算出一個剪切黏度值，最終計算值取這三個值的平均數(shù)。

剪切黏度的計算在NVT系綜下進行，自相關(guān)函數(shù)的相關(guān)長度取8ps，總共計算1000萬步。

剪切黏度計算值與實驗值的比較見圖7，其中紅色的曲線為文獻值，藍色的曲線為模擬值的擬合曲線，擬合曲線與文獻值的平均誤差小于5%。擬合得到的黏度-溫度關(guān)聯(lián)式如下：

2.5 熱導(dǎo)率

熱導(dǎo)率同樣采用平衡分子動力學(xué)計算方法，根據(jù)Green-Kubo公式，由模擬體系熱流通量自相關(guān)函數(shù)（Heat Flux Auto-Correlation Function）計算出體系的熱導(dǎo)率，如式(8)所示。式中B為波爾茲曼常數(shù)，為體系的溫度，Ex為熱流分量E的方向分量。熱流通量在、、方向上的每一個分量都可以計算出一個熱導(dǎo)率值，最終計算值取這3個值的平均數(shù)。

熱導(dǎo)率的計算在NVT系綜下進行，自相關(guān)函數(shù)的相關(guān)長度取8ps，共計算1000萬步。

熱導(dǎo)率的計算值與文獻值的對比見圖8，其中紅色實線為Serrano-López[12]總結(jié)的數(shù)據(jù)庫中的推薦值，虛線為其它文獻[13-15]中報道的熱導(dǎo)率值，可見不同文獻報道的熱導(dǎo)率的值差異很大，本文得到的計算值與Serrano-López所總結(jié)的推薦值之間的平均偏差約24%，考慮到本文所采用的經(jīng)典分子動力學(xué)方法，此誤差在可接受范圍內(nèi)。

2.6 比熱容

比熱容通過體系總能量的漲落來計算，如式(9)～(10)所示。恒容比熱容在NVT系綜下進行，恒壓比熱容在NPT系綜下進行。式中(2)=(2)-()2，為體系總能量?？偣灿嬎?00萬步。

比熱容計算值見表4，將不同溫度下的C進行平均之后得到熔鹽的比熱容為1428.141 J/(kg·K)，文獻值為1561 J/(kg·K)，計算值與文獻值相比偏差為8.5%，兩者吻合得較好。

1.3 常規(guī)治療 術(shù)后給予止血、預(yù)防感染、補液等治療，補液量在3000 mL左右，采用頭低腳高位，頭偏向患側(cè)，術(shù)后第3天給予阿托伐他汀鈣片20 mg口服，每日1次。根據(jù)引流情況術(shù)后3 d內(nèi)拔除引流管。

表4 熔鹽的定壓比熱容計算值

3 結(jié) 論

本文利用Buckingham勢函數(shù)對NaNO3-KNO3- NaNO2三元混合熔鹽進行了分子動力學(xué)計算。徑向分布函數(shù)、配位數(shù)與角分布函數(shù)計算結(jié)果表明，隨著溫度的升高，離子間相互作用減弱，離子團簇變得松散，且異號離子之間距離變大，而同號離子間距沒有明顯變化。配位數(shù)與角分布函數(shù)表明混合熔鹽中陰陽離子呈六配位八面體構(gòu)型，且Na-N的八面體構(gòu)型比K-N更為顯著。同時本文對熔鹽密度、剪切黏度、熱導(dǎo)率、比熱容進行了模擬計算，并與文獻值進行了對比。結(jié)果表明，計算值與文獻值均吻合得很好，驗證了本文所采用的勢函數(shù)與計算方法的可靠性，為下一步大規(guī)模計算模擬熔鹽配方提供了研究基礎(chǔ)。

[1] 國家能源局. 太陽能發(fā)展“十三五”規(guī)劃[R]. 2016.

[2] VIGNAROOBAN K, XU X, ARVAY A, et al. Heat transfer fluids for concentrating solar power systems—A review[J]. Applied Energy, 2015, 146: 383-396.

[3] NUNES V M B, QUEIRóS C S, LOUREN?O M J V, et al. Molten salts as engineering fluids—A review[J]. Applied Energy, 2016, 183: 603-611.

[4] LYNDEN-BELL R M, IMPEY R W, KLEIN M L. Investigation of the lattice vibrations of solid NaNO2by means of molecular dynamics calculations[J]. Chemical Physics, 1986, 109(1): 25-33.

[5] PENG Q, DING J, WEI X, et al. First-principles study for thermodynamic properties of solid KNO2system[J]. International Journal of Thermophysics, 2015, 36(10/11): 2833-2844.

[6] JAYARAMAN S, THOMPSON A P, VON LILIENFELD O A, et al. Molecular simulation of the thermal and transport properties of three alkali nitrate salts[J]. Industrial & Engineering Chemistry Research, 2010, 49(2): 559-571.

[7] URAHATA S R M, RIBEIRO M C C. Molecular dynamics simulation of molten LiNO3with flexible and polarizable anions[J]. Physical Chemistry Chemical Physics, 2003, 5(12): 2619-2624.

[8] RIBEIRO M C C. On the Chemla effect in molten alkali nitrates[J]. The Journal of Chemical Physics, 2002, 117(1): 266-276.

[9] JANZ G J, JAMES D W. Raman spectra and ionic interactions in molten nitrates[J]. The Journal of Chemical Physics, 1961, 35(2): 739-745.

[10] HISATSUNE I C, DEVLIN J P, CALIFANO S. Urey-Bradley potential constants in nitrogen dioxide, nitrite and dinitrogen tetroxide[J]. Spectrochimica Acta, 1960, 16: 450-458.

[11] WANG J, SUN Z, LU G, et al. Molecular dynamics simulations of the local structures and transport coefficients of molten alkali chlorides[J]. The Journal of Physical Chemistry B, 2014, 118(34): 10196-10206.

[12] SERRANO-LóPEZ R, FRADERA J, CUESTA-LóPEZ S. Molten salts database for energy applications[J]. Chemical Engineering and Processing: Process Intensification, 2013, 73: 87-102.

[13] JANZ G J, TOMKINS R P T. Physical properties data compilations relevant to energy storage. IV. Molten salts: Data on additional single and multi-component salt systems[R]. U.S. Department of Commerce, 1981.

[14] WU Y, CHEN C, LIU B, et al. Investigation on forced convective heat transfer of molten salts in circular tubes[J]. International Communications in Heat and Mass Transfer, 2012, 39(10): 1550-1555.

[15] OMOTANL T, NAGASHLMA A. Thermal conductivity of molten salts, HTS and the LiNO3-NaNO3system, using a modified transient hot-wire method[J]. J. Chem. Eng. Data, 1984, 29: 1-3.

Molecular dynamics simulation of structure and physical properties of NaNO3-KNO3-NaNO2ternary phase-change molten salts

NI Haiou, SUN Ze, LU Guimin, YU Jianguo

(National engineering research center for comprehensive utilization of salt lake resources, East China University of Science and Technology, Shanghai 200237, China)

Phase-change salt is an important heat transfer and storage material for concentrated solar power plants. This calls for a quantitative understanding of the relationship between the structure and properties of molten salts. We performed molecular dynamics simulations on NaNO3-KNO3-NaNO2ternary molten salts using Buckingham potential with an aim to understand such relationship. The simulations gave structural information such as radical distribution function, coordination number and angular distribution function, and physical properties including density, shear viscosity, thermal conductivity and heat capacity. The results showed that all the properties calculated agree well with the literature data, suggesting the reliability of pair potential and simulation method used in the work.

molecular dynamics simulation; structure; physical property

10.12028/j.issn.2095-4239.2017.0060

TB 34

A

2095-4239（2017）04-669-06

2017-05-03；

2017-05-16。

國家自然科學(xué)基金（U1407126）；青海省應(yīng)用基礎(chǔ)研究（2017- ZJ-727）；青海省重大科技專項（2013-G-A1A-3）。

倪海歐（1988—），男，主要研究方向為相變材料，E-mail：nihaiou@foxmail.com；

孫澤，博士，副教授，研究方向為熔鹽儲能，化工過程模擬，E-mail：zsun@ecust.edu.cn；路貴民，教授，研究方向為熔鹽離子結(jié)構(gòu)，鹽湖資源綜合利用等，E-mail：gmlu@ecust.edu.cn。