低溫固定點(diǎn)氬的分子力場(chǎng)參數(shù)模擬優(yōu)化方法

摘要：為解決分子力場(chǎng)勢(shì)能函數(shù)中分子間相互作用的Lennard-Jones參數(shù)難以準(zhǔn)確求取的問題，以1990年國際溫標(biāo)所定義的低溫固定點(diǎn)工質(zhì)氬為例，提出一種基于勢(shì)能參數(shù)-臨界特性量效規(guī)律的分子力場(chǎng)Lennard-Jones參數(shù)優(yōu)化方法。基于吉布斯系綜蒙特卡羅方法，利用Trappe-ua和COMPASS兩種力場(chǎng)對(duì)氬的氣液相平衡開展了分子模擬研究，得到初步模擬結(jié)果。通過參數(shù)擾動(dòng)法進(jìn)行兩勢(shì)能參數(shù)對(duì)飽和氣液相密度、飽和蒸氣壓和蒸發(fā)焓影響的定性分析，進(jìn)一步通過定量分析構(gòu)建了COMPASS力場(chǎng)勢(shì)能參數(shù)與臨界特性間的函數(shù)關(guān)系，將實(shí)驗(yàn)臨界值代入確定了優(yōu)化的COMPASS力場(chǎng)的Lennard-Jones參數(shù)。利用優(yōu)化的參數(shù)重新編寫力場(chǎng)進(jìn)行模擬驗(yàn)證。結(jié)果表明，對(duì)于COMPASS力場(chǎng)而言，氬的飽和氣液相密度、蒸氣壓及蒸發(fā)焓模擬精度顯著提升，與Refprop 9.1 軟件計(jì)算值的平均絕對(duì)相對(duì)偏差由27%降到了18%，可為低溫區(qū)其他固定點(diǎn)工質(zhì)氣液相平衡的分子力場(chǎng)Lennard-Jones參數(shù)優(yōu)化提供參考。

關(guān)鍵詞：氬；低溫固定點(diǎn)；氣液相平衡；吉布斯系綜蒙特卡羅方法

中圖分類號(hào)：TK121 文獻(xiàn)標(biāo)志碼：A

DOI：10.7652/xjtuxb202505011 文章編號(hào)：0253-987X（2025）05-0107-13

Simulation and Optimization of Molecular Force Field Parameters

for Low-Temperature Fixed Point Argon

JIA Huilun^1，2，3， LIU Siqi^1，2，4， ZHANG Haiyang^1，2，3， GAO Bo^1，2，3

（1. State Key Laboratory of Cryogenic Science and Technology， Technical Institute of Physics and Chemistry， Chinese Academy of Sciences， Beijing 100190， China; 2.TIPC-LNE Joint Laboratory on Cryogenic Metrology Science and Technology， Technical Institute of Physics and Chemistry， Chinese Academy of Sciences， Beijing 100190， China; 3. University of Chinese Academy of Sciences， Beijing 100490， China; 4. School of Metallurgy， Northeastern University， Shenyang 110819， China）

Abstract：To address the challenge of accurately determining the Lennard-Jones parameters for intermolecular interactions in the potential energy function of molecular force fields， argon， a low-temperature fixed-point working substance defined by the 1990 International Temperature Scale， was taken as an example to propose a molecular force field Lennard-Jones parameter optimization method based on the potential energy parameter-critical property quantity efficiency rule. Initially， molecular simulation studies on the vapor-liquid phase equilibrium of argon were conducted using the Trappe-ua and COMPASS force fields based on the Gibbs ensemble Monte Carlo method， yielding preliminary simulation results. Subsequently， a qualitative analysis was carried out to determine the impact of two potential energy parameters on the saturated vapor-liquid phase density， saturated vapor pressure， and evaporation enthalpy with the parameter perturbation method. In addition， a quantitative analysis was also carried out to establish the functional relationship between the COMPASS force field potential energy parameters and critical properties， and the optimized Lennard-Jones parameters of the COMPASS force field were determined by substituting experimental critical values. Finally， the optimized parameters were used to reprogram the force field for simulation verification. As revealed by the results， for the Compass force field， the average absolute relative deviation of the simulated values of saturated vapor-liquid phase density， vapor pressure， and evaporation enthalpy of argon from the values calculated by the Refprop 9.1 software decreased from 27% to 18%， indicating a significant improvement in simulation accuracy. The method proposed in this study can provide a reference for the optimization of Lennard-Jones parameters for vapor-liquid phase equilibrium of other fixed-point working fluids in the low-temperature region.

Keywords：argon; low-temperature fixed point; vapor-liquid phase equilibrium; Gibbs ensemble Monte Carlo method

低溫固定點(diǎn)作為國際溫標(biāo)^［1］的重要組成部分，在熱力學(xué)溫度賦值及傳遞過程中具有重要意義。當(dāng)前針對(duì)低溫固定點(diǎn)的實(shí)驗(yàn)研究^［2-6］存在周期長(zhǎng)、難度大、成本昂貴等問題，且實(shí)驗(yàn)中無法探究工質(zhì)在相平衡過程中的具體相行為，尤其對(duì)于三相點(diǎn)及臨界點(diǎn)等相變頻繁的區(qū)域無法針對(duì)實(shí)驗(yàn)結(jié)果給出具體解釋，這是當(dāng)前低溫固定點(diǎn)復(fù)現(xiàn)面臨的難題。然而，分子模擬技術(shù)^［7］可以從微觀角度探究低溫固定點(diǎn)工質(zhì)在相平衡狀態(tài)時(shí)的分子構(gòu)型，從而構(gòu)建微觀結(jié)構(gòu)與宏觀熱物理性質(zhì)間的橋梁，且伴隨著計(jì)算機(jī)技術(shù)的發(fā)展，分子模擬與算法的深度耦合^［8-9］使得其模擬準(zhǔn)確性越來越高，這為解決上述難題提供了新的思路。

在固定點(diǎn)工質(zhì)三相平衡中氣液相平衡是其重要的組成部分，當(dāng)前吉布斯系綜蒙特卡羅方法（GEMC）^［10-12］在氣液相平衡分子模擬領(lǐng)域具有顯著優(yōu)勢(shì)，此方法可以不考慮氣液相邊界，不需要直接計(jì)算化學(xué)勢(shì)，且模擬計(jì)算簡(jiǎn)單高效，因此廣泛應(yīng)用于工業(yè)流體及制冷工質(zhì)的氣液相平衡模擬研究。Zhang等^［13］使用GEMC方法研究乙烯丙烯酸-乙酸和乙烯醋酸乙酯-乙酸的氣液平衡；Dong等^［14］利用GEMC研究乙烯醋酸酯系統(tǒng)的氣液相平衡；Zhang等^［15］利用GEMC方法對(duì)HFC-161及其混合物HFC-161+HFO-1234yf的氣液平衡性質(zhì)進(jìn)行分子模擬研究等。為了得到準(zhǔn)確的氣液相平衡模擬結(jié)果，高精度的分子力場(chǎng)是必要的，然而目前針對(duì)低溫固定點(diǎn)的分子力場(chǎng)開發(fā)研究鮮見報(bào)道，Yang等^［16］根據(jù)實(shí)驗(yàn)數(shù)據(jù)采用經(jīng)驗(yàn)擬合法構(gòu)建了包含氬、氧、氖、氫的COMPASS力場(chǎng)，但時(shí)間久遠(yuǎn)并受限于當(dāng)時(shí)的實(shí)驗(yàn)條件及擬合方法，該力場(chǎng)吉布斯系綜蒙特卡羅算法計(jì)算結(jié)果的精度有待提高，而Martin等構(gòu)建的Trappe-ua力場(chǎng)^［17］僅包含了氬的力場(chǎng)參數(shù)，因此完整的低溫固定點(diǎn)的高精度力場(chǎng)仍然缺失。作為分子力場(chǎng)函數(shù)的組成部分，分子間范德華相互作用的Lennard-Jones（LJ）參數(shù)對(duì)于高精度力場(chǎng)的構(gòu)建十分重要。以往獲取LJ參數(shù)的方法主要依賴于對(duì)初始參數(shù)的反復(fù)試錯(cuò)調(diào)整以求得最優(yōu)解，計(jì)算量龐大且耗費(fèi)時(shí)間成本^［18-19］。Bourasseau等提出的LJ參數(shù)迭代優(yōu)化程序復(fù)雜，計(jì)算過程涉及大量偏導(dǎo)，計(jì)算難度大^［20］。采用第一性原理計(jì)算推導(dǎo)LJ參數(shù)時(shí)對(duì)分子間弱相互作用難以進(jìn)行準(zhǔn)確預(yù)測(cè)，且面臨計(jì)算量龐大的問題。

基于此，本文以低溫固定點(diǎn)工質(zhì)氬為例，通過揭示勢(shì)能參數(shù)與臨界特性的量效規(guī)律，提出一種簡(jiǎn)單有效的分子力場(chǎng)LJ參數(shù)優(yōu)化方法，并確定了最佳的LJ勢(shì)能函數(shù)形式，可為低溫區(qū)其他固定點(diǎn)工質(zhì)的LJ參數(shù)優(yōu)化提供參考。

1 Ar的初始模型構(gòu)建及初步模擬結(jié)果

1.1 Ar的初始模型構(gòu)建

本文所有模擬計(jì)算均在開源軟件MCCCS Towhee^［21］中完成，所采用的GEMC方法原理如圖1所示，其采用兩個(gè)熱力學(xué)相關(guān)但相互獨(dú)立的盒子代表氣液兩相，模擬時(shí)通過盒內(nèi)粒子運(yùn)動(dòng)、盒子體積漲落及盒間分子交換3類隨機(jī)運(yùn)動(dòng)分別達(dá)到內(nèi)部平衡、壓力平衡及化學(xué)勢(shì)平衡，從而達(dá)到熱力學(xué)平衡，結(jié)合統(tǒng)計(jì)力學(xué)理論最終得到工質(zhì)熱力學(xué)特性。

針對(duì)低溫固定點(diǎn)工質(zhì)Ar，本文利用Trappe-ua聯(lián)合原子力場(chǎng)^［22］和COMPASS全原子力場(chǎng)^［23］對(duì)其開展了三相點(diǎn)至臨界點(diǎn)溫區(qū)的氣液相平衡分子模擬研究。作為單原子分子，其力場(chǎng)勢(shì)能模型只考慮分子間的范德華作用，其中Trappe-ua力場(chǎng)采用Lennard-Jones 12-6作用勢(shì)（LJ 12-6），COMPASS力場(chǎng)采用Lennard-Jones 9-6作用勢(shì)（LJ 9-6），表達(dá)式如下

E₁=∑ni=1∑nj≥i4ε_ijσ_ijr_ij¹²－σ_ijr_ij⁶（1）

E₂=∑ni=1∑nj≥iε_ij2σ_ijr_ij⁹－3σ_ijr_ij⁶（2）

式中：ε_ij表示兩個(gè)分離的原子對(duì)的勢(shì)能阱深度；σ_ij表示相互作用勢(shì)能正好為0時(shí)的兩原子對(duì)間的距離；r_ij表示兩分子之間實(shí)際距離。

本文模擬時(shí)的初始力場(chǎng)參數(shù)如表1所示，來源于文獻(xiàn)［16-17］，k_b為玻爾茲曼常數(shù)。

根據(jù)模擬過程中參數(shù)的設(shè)置條件不同，GMEC方法分為正則系綜NVT-GEMC和等溫等壓系綜NPT-GEMC方法，本文開展的三相點(diǎn)至臨界點(diǎn)溫區(qū)Ar的純質(zhì)氣液相平衡更適用于NVT-GEMC^［24］。原因在于對(duì)于純質(zhì)的氣液相平衡而言，溫度是影響結(jié)果的主導(dǎo)因素，采用NVT-GEMC方法，在模擬過程中保持粒子的總數(shù)N、盒子的總體積V、模擬溫度T為常數(shù)，更適合開展此類研究。NPT-GEMC可以更好地預(yù)測(cè)不同組分之間的相互作用和平衡狀態(tài)，因此在混合物的氣液相平衡模擬中更加適合使用。本文所有模擬初始均設(shè)液相粒子400個(gè)，氣相粒子100個(gè)，根據(jù)盒內(nèi)粒子數(shù)、模擬溫度及壓力等參數(shù)計(jì)算盒子的初始尺寸，以埃為單位。模擬過程中非鍵作用的截止值依據(jù)盒子的尺寸進(jìn)行調(diào)整，且對(duì)超過截止值的電位部分進(jìn)行長(zhǎng)程尾部分析以進(jìn)行修正。本文所有模擬均進(jìn)行15萬步，其中前10萬步用于體系預(yù)平衡階段，后5萬步用于統(tǒng)計(jì)平衡時(shí)的熱力學(xué)性質(zhì)。每一步循環(huán)都包括了500次蒙特卡羅移動(dòng)，如體積漲落、分子交換、盒內(nèi)分子的運(yùn)動(dòng)，聚集體積偏置移動(dòng)，根據(jù)蔣國柱等^［25］的模擬經(jīng)驗(yàn)，粒子各項(xiàng)運(yùn)動(dòng)的概率占比會(huì)隨模擬溫度進(jìn)行調(diào)整，為了便于觀察模擬的具體過程，模擬階段時(shí)共分了10個(gè)數(shù)據(jù)模塊。

1.2 力場(chǎng)參數(shù)優(yōu)化方法

本文力場(chǎng)參數(shù)的整體優(yōu)化路徑如下所示。

（1）基于Trappe-ua力場(chǎng)與COMPASS力場(chǎng)的初始參數(shù)進(jìn)行NVT-GEMC模擬，分別得到初步模擬結(jié)果。

（2）將初步模擬結(jié)果與Refprop軟件計(jì)算值對(duì)比，選擇優(yōu)化空間較大的力場(chǎng)作為本文優(yōu)化的目標(biāo)，采用參數(shù)擾動(dòng)法進(jìn)行優(yōu)化。參數(shù)擾動(dòng)法即保持一個(gè)參數(shù)不變，控制另一個(gè)參數(shù)等間隔的規(guī)律性變化，并針對(duì)每一組參數(shù)（σ，ε/k_b）在三相點(diǎn)至臨界點(diǎn)溫區(qū)進(jìn)行GEMC模擬計(jì)算，根據(jù)Ar的氣液相平衡特性在整個(gè)溫區(qū)范圍內(nèi)隨可調(diào)參數(shù)變化的模擬結(jié)果，找到兩參數(shù)的影響規(guī)律。

（3）尺寸參數(shù)σ和能量參數(shù)ε/k_b的優(yōu)化可并列進(jìn)行，且優(yōu)化思路是相似的，以下以尺寸參數(shù)σ的優(yōu)化為例。

（4）保持能量參數(shù)ε/k_b不變，尺寸參數(shù)σ等間隔變化：σ+nΔσ，其中Δσ為等間隔變化量，n為變化的份數(shù)。

（5）針對(duì)每一組（σ+nΔσ， ε/k_b）進(jìn)行NVT-GEMC模擬，進(jìn)而得到尺寸參數(shù)σ與模擬所得臨界密度ρ_c的關(guān)系。

（6）將實(shí)驗(yàn)臨界密度ρ_c-exp代入步驟（5）所求得的關(guān)系式中，求得其對(duì)應(yīng)的優(yōu)化后的尺寸參數(shù)。能量參數(shù)的優(yōu)化步驟同理。

（7）基于優(yōu)化參數(shù)進(jìn)行NVT-GEMC模擬，將結(jié)果與Refprop軟件計(jì)算值對(duì)比，驗(yàn)證優(yōu)化效果。

1.3 模擬結(jié)果評(píng)價(jià)指標(biāo)

為評(píng)估力場(chǎng)參數(shù)的模擬效果，本文將模擬結(jié)果與Refprop9.1軟件^［26］計(jì)算值進(jìn)行比較，具體評(píng)價(jià)指標(biāo)如下。Refprop9.1^［26］是由美國國家標(biāo)準(zhǔn)與技術(shù)研究院（NIST）開發(fā)的一款用于計(jì)算制冷工質(zhì)等重要工業(yè)流體熱物理性質(zhì)的軟件數(shù)據(jù)庫，在流體熱物性領(lǐng)域已得到廣泛應(yīng)用。

Ar的飽和氣液相密度的平均絕對(duì)相對(duì)偏差M的計(jì)算式如下

M=∑ni=1|ρ_i，GEMC－ρ_i，Refprop|ρ_i，Refpropn×100%（3）

式中：ρ_i，GEMC是第i個(gè)溫度下飽和氣液相密度的模擬值；ρ_i，Refprop是第i個(gè)溫度下飽和氣液相密度Refprop軟件計(jì)算值；n是計(jì)算的溫度點(diǎn)數(shù)，飽和蒸氣壓與蒸發(fā)焓的平均絕對(duì)相對(duì)偏差計(jì)算同理。

Ar的飽和氣液相密度相對(duì)偏差R的計(jì)算式

如下

R=ρ_i，GEMC－ρ_i，Refpropρ_i，Refprop×100%（4）

Ar的飽和氣液相密度絕對(duì)相對(duì)偏差A(yù)的計(jì)算式如下

A=ρ_i，GEMC－ρ_i，Refpropρ_i，Refprop×100%（5）

1.4 Ar的氣液相平衡初步模擬結(jié)果

1.4.1 飽和氣液相密度

表2列出基于初始力場(chǎng)參數(shù)得到的氬飽和氣液相密度GEMC模擬值與Refprop軟件計(jì)算值，其中：ρ_LR、ρ_LT、ρ_LC分別為Refprop、Trappe-ua力場(chǎng)以及COMPASS力場(chǎng)所得飽和液相密度計(jì)算值；ρ_VR、ρ_VT、ρ_VC分別為Refprop、Trappe-ua力場(chǎng)以及COMPASS力場(chǎng)所得飽和氣相密度計(jì)算值。由表2的數(shù)據(jù)計(jì)算模擬值與Refprop軟件計(jì)算值之間的偏差可得，在整個(gè)溫區(qū)范圍內(nèi)，Trappe-ua力場(chǎng)飽和液相與氣相密度的模擬值相對(duì)于Refprop軟件計(jì)算值的平均絕對(duì)相對(duì)偏差分別是0.50%和6.69%，而COMPASS力場(chǎng)模擬結(jié)果的平均絕對(duì)相對(duì)偏差分別為4.43%和26.62%，表明Trappe-ua力場(chǎng)的飽和氣液相密度模擬結(jié)果與Refprop軟件計(jì)算值吻合更好。其次，兩力場(chǎng)在145 K或85 K附近模擬偏差增大，是因?yàn)榻咏R界點(diǎn)或三相點(diǎn)溫度時(shí)，氣液相變頻繁、分子交換及熱運(yùn)動(dòng)劇烈使得系統(tǒng)難以達(dá)到取樣平均進(jìn)而造成偏差增大。氣相密度模擬偏差明顯大于液相，這是由于凝聚相分子間相互作用較強(qiáng)，使得力場(chǎng)可以做出更準(zhǔn)確的預(yù)測(cè)，同時(shí)氣相密度遠(yuǎn)小于液相密度，使得較小的數(shù)值波動(dòng)可以造成較大的相對(duì)偏差。

1.4.2 飽和蒸氣壓

Ar飽和蒸氣壓GEMC模擬值與Refprop軟件計(jì)算值如表3所示，其中p_R、p_T、p_C分別為Refprop、Trappe-ua力場(chǎng)以及COMPASS力場(chǎng)所得飽和蒸氣壓。基于表3結(jié)果計(jì)算模擬值與Refprop軟件計(jì)算值間的偏差，其中Trappe-ua力場(chǎng)飽和蒸氣壓的模擬值相對(duì)于Refprop軟件計(jì)算值的平均絕對(duì)相對(duì)偏差為6.48%，COMPASS力場(chǎng)為20.58%，因此在Ar的飽和蒸氣壓預(yù)測(cè)方面，Trappe-ua力場(chǎng)更具優(yōu)勢(shì)。表3中計(jì)算結(jié)果表明，COMPASS力場(chǎng)飽和蒸氣壓的模擬值均顯著低于Refprop軟件計(jì)算值，表明其力場(chǎng)勢(shì)能參數(shù)存在較大的優(yōu)化空間。

1.4.3 蒸發(fā)焓

Ar蒸發(fā)焓GEMC模擬值與Refprop軟件計(jì)算值如表4所示，其中ΔH_R、ΔH_T、ΔH_C分別為Refprop、Trappe-ua力場(chǎng)以及COMPASS力場(chǎng)所得蒸發(fā)焓；基于表4結(jié)果計(jì)算模擬值與Refprop軟件計(jì)算值間的偏差，其中Trappe-ua力場(chǎng)蒸發(fā)焓的模擬值與Refprop軟件計(jì)算值的平均絕對(duì)相對(duì)偏差為1.69%，COMPASS力場(chǎng)為12.12%，表明在Ar的蒸發(fā)焓特性計(jì)算方面，Trappe-ua力場(chǎng)更具優(yōu)勢(shì)。蒸發(fā)焓與凝聚相性質(zhì)相關(guān)，表4計(jì)算結(jié)果進(jìn)一步說明，COMPASS力場(chǎng)對(duì)Ar的液相分子相互作用表征方面存在不足。

1.4.4 臨界特性

基于上述模擬結(jié)果，通過density scaling法則^［27］、rectilinear diameters法則^［27］以及Clausius-Clapeyron法則^［28］表達(dá)式如下

ρ_l－ρ_v=A1－TT_c^β（6）

ρ_v+ρ_l2=ρ_c+B1－TT_c（7）

lnp=C₁+C₂T（8）

式中：ρ_v、ρ_l為三相點(diǎn)至臨界點(diǎn)溫區(qū)擬合所得的飽和氣液相密度；p為三相點(diǎn)至臨界點(diǎn)溫區(qū)擬合所得的飽和蒸氣壓；T_c為臨界溫度；ρ_c為臨界密度；A、B、C₁、C₂為常數(shù)；β為臨界因子，此處取0.325。

為了從臨界性質(zhì)的預(yù)測(cè)角度對(duì)力場(chǎng)精度進(jìn)行進(jìn)一步比較，對(duì)臨界溫度、臨界密度及臨界壓力p_c進(jìn)行非線性擬合并與實(shí)驗(yàn)值進(jìn)行對(duì)比。Ar的臨界特性實(shí)驗(yàn)值與GEMC模擬結(jié)果對(duì)比見表5。基于其計(jì)算二者之間的相對(duì)偏差可得，Trappe-ua力場(chǎng)臨界溫度、臨界密度、臨界壓力的擬合計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的相對(duì)偏差幾乎均在3%以內(nèi)，而COMPASS力場(chǎng)臨界性質(zhì)擬合計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的相對(duì)偏差在6%以內(nèi)，進(jìn)一步證明了Trappe-ua力場(chǎng)在計(jì)算Ar的臨界熱力學(xué)性質(zhì)方面的精確性。

綜上所述，兩力場(chǎng)均能合理地表征Ar的熱力學(xué)性質(zhì)隨溫度變化的趨勢(shì)，但Trappe-ua力場(chǎng)在模擬Ar的熱力學(xué)性質(zhì)方面的精確性顯著優(yōu)于COMPASS力場(chǎng)，說明LJ 12-6勢(shì)的準(zhǔn)確性優(yōu)于LJ 9-6勢(shì)，可為其他固定點(diǎn)工質(zhì)LJ勢(shì)能函數(shù)的選擇提供參考。COMPASS力場(chǎng)具有較大的優(yōu)化空間，則其力場(chǎng)參數(shù)優(yōu)化效果更加顯著，本文旨在為其他低溫固定點(diǎn)工質(zhì)的LJ參數(shù)優(yōu)化提供參考，因此第2節(jié)中主要基于COMPASS力場(chǎng)的初步模擬結(jié)果進(jìn)行了力場(chǎng)參數(shù)對(duì)氣液相平衡特性影響的探究。

2 Ar的力場(chǎng)參數(shù)對(duì)氣液相平衡特性的影響

本文所模擬固定點(diǎn)工質(zhì)Ar（作為單原子分子）的力場(chǎng)模型，該模型只考慮分子間范德華作用，包括尺寸參數(shù)σ和能量參數(shù)ε/k_b，基于此，采用參數(shù)擾動(dòng)法探究了兩參數(shù)在三相點(diǎn)至臨界點(diǎn)溫區(qū)范圍內(nèi)對(duì)Ar氣液相平衡的影響規(guī)律。根據(jù)本文此前的模擬經(jīng)驗(yàn)及反復(fù)試錯(cuò)，在初始文獻(xiàn)參數(shù)的基礎(chǔ)上針對(duì)尺寸參數(shù)σ和能量參數(shù)ε/k_b分別以0.2 ?和2 K的等間隔調(diào)動(dòng)，具體結(jié)果如下。

2.1 尺寸參數(shù)σ對(duì)Ar氣液相平衡特性的影響

圖2包含Ar在三相點(diǎn)至臨界點(diǎn)溫區(qū)范圍內(nèi)飽和氣液相密度、蒸氣壓及蒸發(fā)焓隨尺寸參數(shù)σ等間隔變化的模擬結(jié)果，根據(jù)不同尺寸參數(shù)σ進(jìn)行模擬計(jì)算所得的臨界值如表6所示，尺寸參數(shù)的變化對(duì)Trappe-ua力場(chǎng)和COMPASS力場(chǎng)的影響規(guī)律是近似的。圖2（a）和圖2（b）表明，隨著尺寸參數(shù)σ的增

大，飽和氣液相密度曲線整體呈向左平移的趨勢(shì)，即整個(gè)溫區(qū)范圍內(nèi)氣液相密度均在減小，臨界密度逐漸減小，但臨界溫度幾乎未發(fā)生變化，圖2（c）和圖2（d）表明，飽和蒸氣壓曲線隨尺寸參數(shù)的增大整體呈下移趨勢(shì)，同一溫度下的蒸氣壓隨尺寸參數(shù)的增大而減小，且臨界壓力同樣隨尺寸參數(shù)的增大而減小。楊智在針對(duì)制冷工質(zhì)R134開展分子力場(chǎng)與狀態(tài)方程耦合的研究中也發(fā)現(xiàn)類似規(guī)律^［30］。圖2（e）和圖2（f）顯示不同尺寸參數(shù)模擬的蒸發(fā)焓與整體平均值的相對(duì)偏差幾乎都在0.5%以內(nèi)，表明尺寸參數(shù)的變化在整個(gè)溫區(qū)范圍內(nèi)對(duì)蒸發(fā)焓幾乎沒有影響。

相對(duì)偏差計(jì)算如下式

Δ_Trappe-ua=Δ_COMPASS=ΔH_σi－ΔH_avgΔH_avg×100%（9）

ΔH_avg=∑5i=1ΔH_σi5（10）

式中：ΔH_σi表示各尺寸參數(shù)對(duì)應(yīng)的蒸發(fā)焓；ΔH_avg表示各尺寸參數(shù)蒸發(fā)焓的平均值；Δ_Trappe-ua、Δ_COMPASS表示各尺寸參數(shù)對(duì)應(yīng)的蒸發(fā)焓與平均值的相對(duì)偏差。

2.2 能量參數(shù)ε/k_b對(duì)Ar氣液相平衡特性的影響

圖3包含Ar在三相點(diǎn)至臨界點(diǎn)溫區(qū)范圍內(nèi)飽和氣液相密度、蒸氣壓及蒸發(fā)焓隨能量參數(shù)ε/k_b等間隔變化時(shí)的模擬結(jié)果，根據(jù)不同能量參數(shù)ε/k_b進(jìn)行模擬計(jì)算得到的臨界值如表7所示，能量參數(shù)ε/k_b的變化對(duì)兩力場(chǎng)的影響規(guī)律同樣是相似的。圖3（a）和圖3（b）表明，隨著能量參數(shù)的增大，飽和氣液相密度曲線整體上移，整個(gè)溫區(qū)范圍內(nèi)氣相密度減小但液相密度增大，臨界密度幾乎保持不變，但臨界溫度逐漸增大，圖3（c）和圖3（d）表明，隨著能量參數(shù)增大，飽和蒸氣壓曲線整體下移，即同一溫度下蒸氣壓隨能量參數(shù)的增大而減小，但臨界壓力值逐漸增大。圖3（e）和圖3（f）表明蒸發(fā)焓曲線隨能量參數(shù)的增大整體上移，即同一溫度下的蒸發(fā)焓逐漸增大。

3 Ar的力場(chǎng)參數(shù)的優(yōu)化及驗(yàn)證

3.1 Ar的力場(chǎng)參數(shù)優(yōu)化

1.4節(jié)中初始模擬結(jié)果已表明，Trappe-ua力場(chǎng)的精確性明顯優(yōu)于COMPASS力場(chǎng)，則后者的優(yōu)化空間較大，因此本節(jié)主要針對(duì)COMPASS力場(chǎng)的勢(shì)能參數(shù)開展優(yōu)化研究。由第2節(jié)中研究規(guī)律可知，尺寸參數(shù)σ主要影響的是臨界密度ρ_c，對(duì)臨界溫度幾乎沒有影響，能量參數(shù)ε/k_b主要影響臨界溫度T_c，對(duì)臨界密度幾乎沒有影響。因此，理論上，可以通過對(duì)初始力場(chǎng)參數(shù)（σ₁， ε₁/k_b）連續(xù)進(jìn)行 （σ₁， ε₁/k_b+Δε），（σ₁+Δσ， ε₁/k_b+Δε）兩次變化從而達(dá)到優(yōu)化力場(chǎng)參數(shù)的效果，如圖4所示。

通過試錯(cuò)法確定最佳變量（Δσ， Δε）的計(jì)算量較大，本文通過表6中的數(shù)據(jù)擬合尺寸參數(shù)σ與臨界密度ρ_c關(guān)系

ρ_c=166.036σ²－1 692.211σ+4 589.487（11）

其中擬合均方根誤差R_MSE=1.15，決定系數(shù)R₂=0.99， 表明擬合結(jié)果良好。將實(shí)驗(yàn)臨界密度535.60 kg·m^－3代入擬合函數(shù)即可得到最佳尺寸參數(shù)σ_exp的值為3.85，如圖5（a）所示；通過表7中的數(shù)據(jù)可知，能量參數(shù)94.64對(duì)應(yīng)的臨界溫度為150.69 K，近似等于實(shí)驗(yàn)臨界溫度值150.687 K，考慮到擬合誤差，本文將94.64作為最佳能量參數(shù)ε_exp/k_b，如圖5（b）所示。

綜上所述，基于實(shí)驗(yàn)臨界值與所得規(guī)律通過內(nèi)插法得到的Ar的COMPASS力場(chǎng)的優(yōu)化力場(chǎng)參數(shù)分別為3.85 ?與94.64 K。

3.2 優(yōu)化結(jié)果驗(yàn)證

基于優(yōu)化的COMPASS力場(chǎng)進(jìn)行了NVT-GEMC分子模擬，模擬條件設(shè)置與1.1節(jié)相同。結(jié)果表明飽和氣液相密度、飽和蒸氣壓和蒸發(fā)焓整體上都得到了明顯優(yōu)化。COMPASS力場(chǎng)優(yōu)化前后Ar的飽和氣液相密度曲線對(duì)比如圖6所示，圖6（a）表明優(yōu)化后的飽和氣液相密度曲線與Refprop軟件計(jì)算值的一致性更好，圖6（b）表明優(yōu)化后氣液相密度的偏差明顯降低。優(yōu)化前后氣液相密度的模擬數(shù)值如表8所示，其中ρ_LCO、ρ_VCO分別指優(yōu)化后的COMPASS力場(chǎng)所得的飽和液相密度和飽和氣相密度。由表8結(jié)果計(jì)算模擬值與Refprop軟件值間的偏差可得，液相密度與Refprop軟件計(jì)算值的平均絕對(duì)相對(duì)偏差由4.43%優(yōu)化至0.65%，氣相密度由26.62%優(yōu)化至17.76%，優(yōu)化效果顯著。其中，高溫區(qū)尤其是臨界溫區(qū)的優(yōu)化效果最佳，是因?yàn)閮?yōu)化時(shí)采取實(shí)驗(yàn)臨界溫度與實(shí)驗(yàn)臨界密度作為目標(biāo)值，因此越接近臨界溫度區(qū)域優(yōu)化效果越好。在三相點(diǎn)溫度附近氣相密度的相對(duì)偏差較大是因?yàn)榇藭r(shí)氣相密度本身較小，因此即使較小的絕對(duì)偏差也會(huì)導(dǎo)致相對(duì)偏差很大，這也進(jìn)一步導(dǎo)致了從相對(duì)偏差角度來看，氣相密度優(yōu)化效果不明顯。

COMPASS力場(chǎng)優(yōu)化前后Ar的飽和蒸氣壓與蒸發(fā)焓曲線對(duì)比如圖7（a）和圖7（b）所示，模擬值如表9所示，其中p_CO、ΔH_CO分別指優(yōu)化后的COMPASS力場(chǎng)所得的飽和蒸氣壓和蒸發(fā)焓。其表明優(yōu)化后的飽和蒸氣壓與蒸發(fā)焓曲線與Refprop軟件計(jì)算值一致性更好。優(yōu)化前后與Refprop軟件計(jì)算值絕對(duì)相對(duì)偏差對(duì)比如圖8（a）和圖8（b）所示，其表明優(yōu)化后的偏差明顯降低，其中飽和蒸氣壓與Refprop軟件計(jì)算值的平均絕對(duì)相對(duì)偏差從20.58% 優(yōu)化至16.84%，蒸發(fā)焓由12.12%優(yōu)化至6.54%，且高溫區(qū)優(yōu)化效果最好，原因在于Ar的飽和蒸氣壓在三相點(diǎn)溫度附近同樣較小，此時(shí)較小的壓力變動(dòng)就會(huì)造成較大的相對(duì)偏差，使得蒸氣壓整體平均絕對(duì)相對(duì)偏差的優(yōu)化效果不明顯。

上述研究結(jié)果表明，利用本文提出的力場(chǎng)參數(shù)優(yōu)化方法對(duì)COMPASS力場(chǎng)參數(shù)進(jìn)行優(yōu)化，各項(xiàng)熱力學(xué)性質(zhì)的NVT-GEMC模擬結(jié)果與Refprop軟件計(jì)算值的平均絕對(duì)相對(duì)偏差從優(yōu)化前的27%降至優(yōu)化后的18%，相對(duì)精度提升了33%，表明了本方法的有效性。此外，與傳統(tǒng)LJ參數(shù)優(yōu)化方法相比，本文提出的參數(shù)擾動(dòng)法在優(yōu)化LJ參數(shù)時(shí)計(jì)算量小，節(jié)省時(shí)間成本。

以上探究規(guī)律表明，無論是Trappe-ua力場(chǎng)還是COMPASS力場(chǎng)，當(dāng)只考慮分子間范德華作用時(shí)，尺寸參數(shù)σ和能量參數(shù)ε/k_b的影響規(guī)律是一致的，由此可將本文所探究的規(guī)律推廣到其他任意分子的LJ參數(shù)優(yōu)化，一定程度上拓寬了本文所研究規(guī)律的應(yīng)用范圍。

4 結(jié) 論

本文以1990年國際溫標(biāo)所定義的低溫固定點(diǎn)工質(zhì)Ar為例，明晰了勢(shì)能參數(shù)對(duì)臨界特性的影響規(guī)律，從而提出了一種簡(jiǎn)單有效的分子力場(chǎng)LJ參數(shù)優(yōu)化方法，并確定了最佳的LJ勢(shì)能函數(shù)形式，可為低溫區(qū)其他固定點(diǎn)工質(zhì)的LJ參數(shù)優(yōu)化提供參考，具體結(jié)論如下。

（1）Trappe-ua聯(lián)合原子力場(chǎng)和COMPASS全原子力場(chǎng)的初始文獻(xiàn)力場(chǎng)參數(shù)均能正確預(yù)測(cè)Ar在三相點(diǎn)至臨界點(diǎn)溫區(qū)熱力學(xué)性質(zhì)隨溫度變化的總體趨勢(shì)，Trappe-ua力場(chǎng)各項(xiàng)模擬結(jié)果與Refprop軟件計(jì)算值的平均絕對(duì)相對(duì)偏差均在7%以內(nèi)，相比于COMPASS力場(chǎng)的27%具有明顯優(yōu)勢(shì)，表明LJ 12-6勢(shì)的準(zhǔn)確性優(yōu)于LJ 9-6勢(shì)，可為其他低溫固定點(diǎn)工質(zhì)LJ作用勢(shì)的選取提供參考。

（2）在三相點(diǎn)至臨界點(diǎn)溫區(qū)范圍內(nèi)，Ar的同一種熱力學(xué)特性隨力場(chǎng)勢(shì)能參數(shù)的變化呈相同變化趨勢(shì)。其中：尺寸參數(shù)σ主要影響臨界密度ρ_c與臨界壓力p_c，隨著尺寸參數(shù)σ的增大，臨界密度ρ_c和臨界壓力p_c降低，而蒸發(fā)焓ΔH和臨界溫度值T_c幾乎沒有變化；能量參數(shù)ε/k_b主要影響臨界溫度T_c、臨界壓力p_c和蒸發(fā)焓ΔH，隨著能量參數(shù)ε/k_b的增大，臨界溫度T_c、臨界壓力p_c和蒸發(fā)焓ΔH增大，而臨界密度值ρ_c幾乎沒有變化。

（3）以實(shí)驗(yàn)臨界溫度和實(shí)驗(yàn)臨界密度為優(yōu)化目標(biāo)，基于上述規(guī)律對(duì)COMPASS力場(chǎng)參數(shù)進(jìn)行了優(yōu)化，結(jié)果表明：優(yōu)化后的力場(chǎng)參數(shù)使得Ar的液相密度的平均絕對(duì)相對(duì)偏差由4.43%降低至0.65%，氣相密度的平均絕對(duì)相對(duì)偏差從26.62%降低至17.76%；飽和蒸氣壓平均絕對(duì)相對(duì)偏差由20.58%降低至16.84%，蒸發(fā)焓平均絕對(duì)相對(duì)偏差由12.12% 降低至6.54%，證明了該方法的有效性。

本文所提出的分子力場(chǎng)優(yōu)化方法可簡(jiǎn)單有效地對(duì)LJ參數(shù)進(jìn)行優(yōu)化，可為后續(xù)低溫區(qū)其他固定點(diǎn)工質(zhì)的分子力場(chǎng)LJ參數(shù)優(yōu)化提供參考，進(jìn)而獲得準(zhǔn)確的氣液相平衡模擬數(shù)據(jù)。同時(shí)也為低溫固定點(diǎn)的氣、液、固三相平衡的分子模擬研究奠定基礎(chǔ)，進(jìn)一步為修正固定點(diǎn)熱力學(xué)溫度及構(gòu)建寬溫區(qū)范圍內(nèi)更加準(zhǔn)確的R-T內(nèi)插方程做出貢獻(xiàn)，有助于實(shí)現(xiàn)熱力學(xué)溫度的精確賦值及有效傳遞，更深層次地為深空探測(cè)、量子計(jì)算等眾多低溫領(lǐng)域的前沿科學(xué)研究提供參考。

參考文獻(xiàn)：

［1］PRESTON-THOMAS H. The international temperature scale of 1990 （ITS-90） ［J］. Metrologia， 1990， 27（1）： 3-10.

［2］PAN Changzhao， SPARASCI F， ZHANG Haiyang， et al. Acoustic measurement of the triple point of neon T_Ne and thermodynamic calibration of a transfer standard for accurate cryogenic thermometry ［J］. Metrologia， 2021， 58（4）： 045006.

［3］SUN Jianping， TSAI S F， LI Ting. Comparison report between NIM and CMS of the ITS-90 realization at the triple point of argon ［J］. Metrologia， 2023， 60（1A）： 03004.

［4］LOPARDO G， DEMATTEIS R， STEUR P P M. Realisation of the triple-point of argon： comparison between two devices ［J］. Physica Scripta， 2024， 99（7）： 075912.

［5］ROURKE P M C. Thermodynamic temperature of the triple point of xenon measured by refractive index gas thermometry ［J］. Metrologia， 2020， 57（2）： 024001.

［6］STEUR P P M， ROURKE P M C， GIRAUDI D. Comparison of xenon triple point realizations ［J］. Metrologia， 2019， 56（1）： 015008.

［7］CICCOTTI G， DELLAGO C， FERRARIO M， et al. Molecular simulations： past， present， and future （a topical issue in EPJB） ［J］. The European Physical Journal： B， 2022， 95（1）： 3.

［8］DOERR S， MAJEWSKI M， PéREZ A， et al. TorchMD： a deep learning framework for molecular simulations ［J］. Journal of Chemical Theory and Computation， 2021， 17（4）： 2355-2363.

［9］LIU Jie， ZHANG Tao， SUN Shuyu. Review of deep learning algorithms in molecular simulations and perspective applications on petroleum engineering ［J］. Geoscience Frontiers， 2024， 15（2）： 101735.

［10］ZHANG Liwen， YANG Yuhong， YIN Kun， et al. A review of GEMC method and its improved algorithms ［J］. Acta Geochimica， 2023， 42（3）： 409-434.

［11］SOROUSH BARHAGHI M， POTOFF J J. Prediction of phase equilibria and Gibbs free energies of transfer using molecular exchange Monte Carlo in the Gibbs ensemble ［J］. Fluid Phase Equilibria， 2019， 486： 106-118.

［12］涂茂萍， 張丹， 袁洋， 等. 靜止顆粒群熱輻射吸收性能的蒙特卡羅法二維數(shù)值研究 ［J］. 西安交通大學(xué)學(xué)報(bào)， 2024， 58（5）： 167-178.

TU Maoping， ZHANG Dan， YUAN Yang， et al. Thermal radiation absorption properties of two-dimensional stationary particle groups based on the Monte-Carlo method ［J］. Journal of Xi’an Jiaotong University， 2024， 58（5）： 167-178.

［13］ZHANG Minhua， WANG Yuankang， LIU Shuhan， et al. A study on the vapor-liquid equilibria of vinyl acrylate-acetic acid and vinyl acetate-acetic acid with GEMC method ［J］. Separation and Purification Technology， 2024， 330 （Part C）： 125509.

［14］DONG Xiuqin， GUAN Xiaoxiao， JIANG Yuan， et al. Extension of the TraPPE-UA force field to the simulation of vapor-liquid phase equilibria of vinyl acetate system ［J］. Journal of Molecular Liquids， 2015， 209： 520-525.

［15］ZHANG Nan， HU Peng， CHEN Longxiang， et al. Molecular modeling of vapor-liquid equilibrium properties of HFC-161 and its mixture HFC-161+HFO-1234yf ［J］. Journal of Molecular Liquids， 2020， 306： 112896.

［16］YANG Jie， REN Yi， TIAN Anmin， et al. COMPASS force field for 14 inorganic molecules， He， Ne， Ar， Kr， Xe， H₂， O₂， N₂， NO， CO， CO₂， NO₂， CS₂， and SO₂， in liquid phases ［J］. The Journal of Physical Chemistry： B， 2000， 104（20）： 4951-4957.

［17］MARTIN M G， SIEPMANN J I. Calculating Gibbs free energies of transfer from Gibbs ensemble Monte Carlo simulations ［J］. Theoretical Chemistry Accounts， 1998， 99（5）： 347-350.

［18］RAABE G， MAGINN E J. A force field for 3， 3， 3-fluoro-1-propenes， including HFO-1234yf ［J］. The Journal of Physical Chemistry： B， 2010， 114（31）： 10133-10142.

［19］YANG Zhi， GONG Maoqiong， DONG Xueqiang， et al. Molecular modeling and simulation of vapor-liquid equilibrium of the refrigerant R152a and its mixture R152a+R32 ［J］. Fluid Phase Equilibria， 2015， 394： 93-100.

［20］BOURASSEAU E， HABOUDOU M， BOUTIN A， et al. New optimization method for intermolecular potentials： optimization of a new anisotropic united atoms potential for olefins： prediction of equilibrium properties ［J］. The Journal of Chemical Physics， 2003， 118（7）： 3020-3034.

［21］CUMMINGS P T， GILMER J B. Open-source molecular modeling software in chemical engineering ［J］. Current Opinion in Chemical Engineering， 2019， 23： 99-105.

［22］JARADAT A， AL-SALMAN R， OBEIDAT A. Molecular dynamics simulation of vapor-liquid equilibrium in 1-alkanol unary systems： a study of surface tension， density， and vapor pressure of TraPPE-UA force field ［J］. Fluid Phase Equilibria， 2024， 577： 113967.

［23］SAVIN A V， MAZO M A. The COMPASS force field： validation for carbon nanoribbons ［J］. Physica： E Low-Dimensional Systems and Nanostructures， 2020， 118： 113937.

［24］SOBECKI N， NIETO-DRAGHI C， DI LELLA A， et al. Phase behavior of hydrocarbons in nano-pores ［J］. Fluid Phase Equilibria， 2019， 497： 104-121.

［25］蔣國柱， 薛榮書， 魏順安. 吉布斯系綜Monte Carlo技術(shù)模擬L-J流體汽液相平衡 ［J］. 重慶大學(xué)學(xué)報(bào)（自然科學(xué)版）， 2003， 26（8）： 64-67.

JIANG Guozhu， XUE Rongshu， WEI Shunan. Computer simulation for liquid-vapour phase equilibrium of L-J fluid by the GEMC method ［J］. Journal of Chongqing University（Natural Science Edition）， 2003， 26（8）： 64-67.

［26］HUBER M L， LEMMON E W， BELL I H， et al. The NIST REFPROP database for highly accurate properties of industrially important fluids ［J］. Industrial amp; Engineering Chemistry Research， 2022， 61（42）： 15449-15472.

［27］YANG Zhiqiang， VALTZ A， COQUELET C， et al. Critical properties and vapor-liquid equilibrium of two near-azeotropic mixtures containing HFOs ［J］. International Journal of Refrigeration， 2022， 138： 133-147.

［28］OHE S. A prediction method of vapor pressures for pure substances ［J］. Fluid Phase Equilibria， 2022， 561： 113474.

［29］GILGEN R， KLEINRAHM R， WAGNERW. Measurement and correlation of the （pressure， density， temperature） relation of argon： I The homogeneous gas and liquid regions in the temperature range from 90 K to 340 K at pressures up to 12 MPa ［J］. The Journal of Chemical Thermodynamics， 1994， 26（4）： 383-398.

［30］楊智. 低GWP工質(zhì)分子力場(chǎng)的開發(fā)及其熱物性的分子模擬 ［D］. 北京： 中國科學(xué)院大學(xué)， 2015： 127-145.

（編輯 武紅江）