低溫固定點氬的分子力場參數模擬優化方法

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

關鍵詞：氬；低溫固定點；氣液相平衡；吉布斯系綜蒙特卡羅方法

中圖分類號：TK121 文獻標志碼：A

DOI：10.7652/xjtuxb202505011 文章編號：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

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

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

基于此，本文以低溫固定點工質氬為例，通過揭示勢能參數與臨界特性的量效規律，提出一種簡單有效的分子力場LJ參數優化方法，并確定了最佳的LJ勢能函數形式，可為低溫區其他固定點工質的LJ參數優化提供參考。

1 Ar的初始模型構建及初步模擬結果

1.1 Ar的初始模型構建

本文所有模擬計算均在開源軟件MCCCS Towhee^［21］中完成，所采用的GEMC方法原理如圖1所示，其采用兩個熱力學相關但相互獨立的盒子代表氣液兩相，模擬時通過盒內粒子運動、盒子體積漲落及盒間分子交換3類隨機運動分別達到內部平衡、壓力平衡及化學勢平衡，從而達到熱力學平衡，結合統計力學理論最終得到工質熱力學特性。

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

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

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

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

本文模擬時的初始力場參數如表1所示，來源于文獻［16-17］，k_b為玻爾茲曼常數。

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

1.2 力場參數優化方法

本文力場參數的整體優化路徑如下所示。

（1）基于Trappe-ua力場與COMPASS力場的初始參數進行NVT-GEMC模擬，分別得到初步模擬結果。

（2）將初步模擬結果與Refprop軟件計算值對比，選擇優化空間較大的力場作為本文優化的目標，采用參數擾動法進行優化。參數擾動法即保持一個參數不變，控制另一個參數等間隔的規律性變化，并針對每一組參數（σ，ε/k_b）在三相點至臨界點溫區進行GEMC模擬計算，根據Ar的氣液相平衡特性在整個溫區范圍內隨可調參數變化的模擬結果，找到兩參數的影響規律。

（3）尺寸參數σ和能量參數ε/k_b的優化可并列進行，且優化思路是相似的，以下以尺寸參數σ的優化為例。

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

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

（6）將實驗臨界密度ρ_c-exp代入步驟（5）所求得的關系式中，求得其對應的優化后的尺寸參數。能量參數的優化步驟同理。

（7）基于優化參數進行NVT-GEMC模擬，將結果與Refprop軟件計算值對比，驗證優化效果。

1.3 模擬結果評價指標

為評估力場參數的模擬效果，本文將模擬結果與Refprop9.1軟件^［26］計算值進行比較，具體評價指標如下。Refprop9.1^［26］是由美國國家標準與技術研究院（NIST）開發的一款用于計算制冷工質等重要工業流體熱物理性質的軟件數據庫，在流體熱物性領域已得到廣泛應用。

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

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

式中：ρ_i，GEMC是第i個溫度下飽和氣液相密度的模擬值；ρ_i，Refprop是第i個溫度下飽和氣液相密度Refprop軟件計算值；n是計算的溫度點數，飽和蒸氣壓與蒸發焓的平均絕對相對偏差計算同理。

Ar的飽和氣液相密度相對偏差R的計算式

如下

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

Ar的飽和氣液相密度絕對相對偏差A的計算式如下

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

1.4 Ar的氣液相平衡初步模擬結果

1.4.1 飽和氣液相密度

表2列出基于初始力場參數得到的氬飽和氣液相密度GEMC模擬值與Refprop軟件計算值，其中：ρ_LR、ρ_LT、ρ_LC分別為Refprop、Trappe-ua力場以及COMPASS力場所得飽和液相密度計算值；ρ_VR、ρ_VT、ρ_VC分別為Refprop、Trappe-ua力場以及COMPASS力場所得飽和氣相密度計算值。由表2的數據計算模擬值與Refprop軟件計算值之間的偏差可得，在整個溫區范圍內，Trappe-ua力場飽和液相與氣相密度的模擬值相對于Refprop軟件計算值的平均絕對相對偏差分別是0.50%和6.69%，而COMPASS力場模擬結果的平均絕對相對偏差分別為4.43%和26.62%，表明Trappe-ua力場的飽和氣液相密度模擬結果與Refprop軟件計算值吻合更好。其次，兩力場在145 K或85 K附近模擬偏差增大，是因為接近臨界點或三相點溫度時，氣液相變頻繁、分子交換及熱運動劇烈使得系統難以達到取樣平均進而造成偏差增大。氣相密度模擬偏差明顯大于液相，這是由于凝聚相分子間相互作用較強，使得力場可以做出更準確的預測，同時氣相密度遠小于液相密度，使得較小的數值波動可以造成較大的相對偏差。

1.4.2 飽和蒸氣壓

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

1.4.3 蒸發焓

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

1.4.4 臨界特性

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

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

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

lnp=C₁+C₂T（8）

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

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

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

2 Ar的力場參數對氣液相平衡特性的影響

本文所模擬固定點工質Ar（作為單原子分子）的力場模型，該模型只考慮分子間范德華作用，包括尺寸參數σ和能量參數ε/k_b，基于此，采用參數擾動法探究了兩參數在三相點至臨界點溫區范圍內對Ar氣液相平衡的影響規律。根據本文此前的模擬經驗及反復試錯，在初始文獻參數的基礎上針對尺寸參數σ和能量參數ε/k_b分別以0.2 ?和2 K的等間隔調動，具體結果如下。

2.1 尺寸參數σ對Ar氣液相平衡特性的影響

圖2包含Ar在三相點至臨界點溫區范圍內飽和氣液相密度、蒸氣壓及蒸發焓隨尺寸參數σ等間隔變化的模擬結果，根據不同尺寸參數σ進行模擬計算所得的臨界值如表6所示，尺寸參數的變化對Trappe-ua力場和COMPASS力場的影響規律是近似的。圖2（a）和圖2（b）表明，隨著尺寸參數σ的增

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

相對偏差計算如下式

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

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

式中：ΔH_σi表示各尺寸參數對應的蒸發焓；ΔH_avg表示各尺寸參數蒸發焓的平均值；Δ_Trappe-ua、Δ_COMPASS表示各尺寸參數對應的蒸發焓與平均值的相對偏差。

2.2 能量參數ε/k_b對Ar氣液相平衡特性的影響

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

3 Ar的力場參數的優化及驗證

3.1 Ar的力場參數優化

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

通過試錯法確定最佳變量（Δσ， Δε）的計算量較大，本文通過表6中的數據擬合尺寸參數σ與臨界密度ρ_c關系

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

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

綜上所述，基于實驗臨界值與所得規律通過內插法得到的Ar的COMPASS力場的優化力場參數分別為3.85 ?與94.64 K。

3.2 優化結果驗證

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

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

上述研究結果表明，利用本文提出的力場參數優化方法對COMPASS力場參數進行優化，各項熱力學性質的NVT-GEMC模擬結果與Refprop軟件計算值的平均絕對相對偏差從優化前的27%降至優化后的18%，相對精度提升了33%，表明了本方法的有效性。此外，與傳統LJ參數優化方法相比，本文提出的參數擾動法在優化LJ參數時計算量小，節省時間成本。

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

4 結 論

本文以1990年國際溫標所定義的低溫固定點工質Ar為例，明晰了勢能參數對臨界特性的影響規律，從而提出了一種簡單有效的分子力場LJ參數優化方法，并確定了最佳的LJ勢能函數形式，可為低溫區其他固定點工質的LJ參數優化提供參考，具體結論如下。

（1）Trappe-ua聯合原子力場和COMPASS全原子力場的初始文獻力場參數均能正確預測Ar在三相點至臨界點溫區熱力學性質隨溫度變化的總體趨勢，Trappe-ua力場各項模擬結果與Refprop軟件計算值的平均絕對相對偏差均在7%以內，相比于COMPASS力場的27%具有明顯優勢，表明LJ 12-6勢的準確性優于LJ 9-6勢，可為其他低溫固定點工質LJ作用勢的選取提供參考。

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

（3）以實驗臨界溫度和實驗臨界密度為優化目標，基于上述規律對COMPASS力場參數進行了優化，結果表明：優化后的力場參數使得Ar的液相密度的平均絕對相對偏差由4.43%降低至0.65%，氣相密度的平均絕對相對偏差從26.62%降低至17.76%；飽和蒸氣壓平均絕對相對偏差由20.58%降低至16.84%，蒸發焓平均絕對相對偏差由12.12% 降低至6.54%，證明了該方法的有效性。

本文所提出的分子力場優化方法可簡單有效地對LJ參數進行優化，可為后續低溫區其他固定點工質的分子力場LJ參數優化提供參考，進而獲得準確的氣液相平衡模擬數據。同時也為低溫固定點的氣、液、固三相平衡的分子模擬研究奠定基礎，進一步為修正固定點熱力學溫度及構建寬溫區范圍內更加準確的R-T內插方程做出貢獻，有助于實現熱力學溫度的精確賦值及有效傳遞，更深層次地為深空探測、量子計算等眾多低溫領域的前沿科學研究提供參考。

參考文獻：

［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］涂茂萍， 張丹， 袁洋， 等. 靜止顆粒群熱輻射吸收性能的蒙特卡羅法二維數值研究 ［J］. 西安交通大學學報， 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技術模擬L-J流體汽液相平衡 ［J］. 重慶大學學報（自然科學版）， 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工質分子力場的開發及其熱物性的分子模擬 ［D］. 北京： 中國科學院大學， 2015： 127-145.

（編輯 武紅江）