韓 勇, 郭向利, 龍新平

(1. 中國工程物理研究院化工材料研究所, 四川 綿陽 621999; 2. 中國工程物理研究院, 四川 綿陽 621999)

1 引 言

炸藥爆轟時,爆轟產物處于高溫高壓狀態(tài),爆轟產物壓力由數(shù)十吉帕逐漸衰減至兆帕量級,溫度也由數(shù)千開爾文逐漸衰減至1000K甚至更低,建立能夠準確描述爆轟產物氣體高溫、高壓的狀態(tài)方程對有效表征爆轟產物狀態(tài)變化過程具有十分重要的意義。

目前,實現(xiàn)爆轟產物氣體所處的高溫、高壓試驗條件十分困難,關于高溫、高壓氣體狀態(tài)的實驗數(shù)據(jù)很少。分子動力學(MD)/蒙特卡洛方法(MC)是研究爆轟產物氣體高溫高壓熱力學性質的有力工具,它們通過有限的分子數(shù)目,結合邊界的約束處理,能夠準確計算氣體的宏觀熱力學性質。然而MD方法和MC方法計算耗時,且無法直接應用。解析形式的狀態(tài)方程則可以直接應用于實際問題,形成子程序直接納入相關程序中應用。在描述爆轟產物氣體方面,已有BKW[1]、VLW[2]、LJD[3]、JCZ[4]等解析狀態(tài)方程形式。其中,VLW狀態(tài)方程為我國吳雄教授借鑒相似理論,認為各階維里系數(shù)在高溫下是相似的,高階維里系數(shù)可以通過二階維里系數(shù)求得,進而將維里物態(tài)方程以一種簡化形式寫出,形成了VLW爆轟產物狀態(tài)方程[5-6],其成功應用于炸藥爆轟性能的理論計算。與BKW狀態(tài)方程相比,VLW狀態(tài)方程在描述爆轟環(huán)境下氣體高溫高壓熱力學狀態(tài)時有了一定改善。然而VLW狀態(tài)方程形式仍然存在一些不足,VLW狀態(tài)方程不能有效描述爆轟產物氣體組分的高溫高壓熱力學狀態(tài),從而使其所采用的爆轟產物勢參數(shù)缺乏必要的物理基礎支撐。作者[7]前期應用VLW狀態(tài)方程計算H2O沖擊Hugoniot曲線的結果表明,即使在優(yōu)化后的勢參數(shù)條件下,VLW狀態(tài)方程仍然不能在寬的壓力范圍有效描述H2O的高溫高壓狀態(tài)。

針對VLW狀態(tài)方程存在的上述不足,本研究提出了一種描述爆轟環(huán)境下高溫、高壓氣體的對比態(tài)維里型狀態(tài)方程VHL(Viral-Han-Long),該狀態(tài)方程基于LJ勢能函數(shù),通過對炸藥爆轟產物中的重要氣體組分CO2高溫高壓熱力學狀態(tài)的描述,表明對比態(tài)VHL狀態(tài)方程能夠有效地描述爆轟環(huán)境下CO2氣體的壓力、體積和溫度(pVT)熱力學關系。

2 對比態(tài)VHL狀態(tài)方程

理論上,任何氣體的狀態(tài)方程,都可以用維里形式描述。然而,實踐中隨著維里系數(shù)階數(shù)的提高,計算的復雜性迅速增大。Barker等[8]基于LJ勢能函數(shù),針對無量綱第三階、第四階、第五階系數(shù),進行了精確的理論計算,其以表格形式列出了特定溫度時的系數(shù)值,不適宜實際使用。在前期研究中,作者等[9]提出用一種簡化維里型狀態(tài)方程形式描述高溫甲烷氣體的熱力學pVT狀態(tài)。為有效描述爆轟環(huán)境下其它氣體產物的pVT熱力學關系,本研究提出了一種基于對比態(tài)原理的維里型狀態(tài)方程形式VHL,表達式見(1)。

(1)

C*=c1T*c2+c3T*c4+c5T*c6+c7T*c8+c9T*c10

D*=d1T*d2+d3T*d4+d5T*d6+d7T*d8+d9T*d10

E*=e1T*e2+e3T*e4

b0=2πNσ3/3

式中,p為氣體實際壓力,GPa;V為摩爾體積,cm3·mol-1;T為溫度,K;R為摩爾氣體常數(shù)8.3145,J·mol-1·K-1;N為阿佛加德羅常數(shù)6.022045×1023mol-1,第二階無量綱維里系數(shù)B*采用變步長辛卜生求積法近似計算獲得[10];C*、D*分別為第三階和第四階無量綱維里系數(shù),c1～c10、d1～d10為常數(shù),通過該表達式所得第三階和第四階無量綱維里系數(shù)值與理論值十分吻合; 對于第五階以上的無量綱維里系數(shù),則采用組合函數(shù)表示,其系數(shù)值e1～e4、a、b、f由無極性分子CH4的熱力學數(shù)據(jù)[11-12](溫度1000 K以上的112組pVT數(shù)據(jù))確定; 各階維里系數(shù)的擬合常數(shù)獲得方式見文獻[9]所述,具體數(shù)值見表1。T*為無量綱溫度,bCH4為67.21,w為對比態(tài)參量,VHL狀態(tài)方程通過參量w實現(xiàn)其它氣體與甲烷狀態(tài)方程的對比。ε、σ為LJ勢參數(shù)。

表1VHL狀態(tài)方程中各階維里系數(shù)擬合常數(shù)值

Table1Fitting constant values of each order viral coefficient in VHL Equation of State(EOS)

coefficientvaluecoefficientvaluecoefficientvaluec1-0.96665 d1-2.31154 e12.15701c2-4.51039 d2-7.49608 e2-2.52736c3-2.58733 d3-2.56864 e30.27080c4-0.85487 d4-1.14947 e4-1.00190c52.08033 d52.09534 a-0.00243c6-0.51631 d6-0.84757 b0.00018c72.02825 d72.81486 f3.60800c8-2.15543 d8-4.44480c9-0.12489 d9-0.30669c10-7.65161 d10-11.36945

3 CO2氣體高溫高壓熱力學狀態(tài)的計算

由于液態(tài)CO2樣品制備條件苛刻、難度大,其高溫、高壓的基礎實驗數(shù)據(jù)較少。1990年初,Schott等[13]才發(fā)表了用化爆技術在5～30 GPa區(qū)域獲得的一組試驗點。Nellis等[14]用二級輕氣炮技術在25～70 GPa區(qū)域獲得了一組試驗點。劉福生等[15]利用二級輕氣炮作沖擊加載手段,獲得了CO2在20～60 GPa區(qū)域六個Hugoniot數(shù)據(jù)點。實驗數(shù)據(jù)點偏少,且缺乏溫度的直接測量數(shù)據(jù)。因此,本研究以Belonoshko等[12]的分子動力學計算數(shù)據(jù)作為CO2的基礎數(shù)據(jù),溫度范圍為718～4978 K,壓力范圍為0.5116～111.078 GPa,基于VHL狀態(tài)方程,應用復形調優(yōu)法[16]優(yōu)化了CO2的LJ勢參數(shù),勢參數(shù)值與文獻值比較見表2所示。采用VHL、VLW狀態(tài)方程計算結果和分子動力學計算值比較見表3所示,不同溫度或壓力下CO2體積計算偏差如圖1、圖2所示。由表3、圖1和圖2可得,采用VHL狀態(tài)方程計算得CO2體積平均絕對偏差為0.971%,最大偏差為4.04%, VHL狀態(tài)方程計算所得體積偏差與壓力和溫度參量無明顯的相關性。VLW狀態(tài)方程計算所得體積偏差則與溫度具有明顯的相關性,在較低溫度下,計算所得體積偏差較大，最大偏差87.149%,隨著溫度的升高,計算體積偏差逐漸減小,但仍普遍高于VHL狀態(tài)方程計算結果,采用VLW狀態(tài)方程計算所得平均絕對偏差20.2%。

表2本研究所采用CO2勢參數(shù)值與文獻值比較

Table2Comparison of the potential parameter values of CO2used in this paper and literature ones

(ε/k)/Kb0/mL·mol-1reference247.063.37[17]205.085.05[5]181.867.30thispaper

圖1不同壓力下VHL狀態(tài)方程、VLW狀態(tài)方程預測CO2的體積誤差百分比

Fig.1Error percentage of volume predicted by VHL EOS and VLW EOS at different pressure

圖2不同溫度下VHL狀態(tài)方程、VLW狀態(tài)方程預測CO2的體積誤差百分比

Fig.2Error percentage of volume predicted by VHL EOS and VLW EOS at different temperature

為驗證對比態(tài)VHL狀態(tài)方程在更低壓力下的有效性,同時與VLW狀態(tài)方程計算結果比較,本研究引

用NIST數(shù)據(jù)庫的數(shù)據(jù)[18],對CO2在1000 K,20～800 MPa的熱力學狀態(tài)進行了計算,結果見表4所示。采用對比態(tài)VHL狀態(tài)方程計算得CO2體積絕對平均偏差為0.698%,采用VLW狀態(tài)方程計算所得體積絕對平均偏差為11.988%。在固定溫度1000 K條件下,隨著壓力的增加,VLW計算所得體積偏差逐漸增大,在壓力為800 MPa時,體積偏差最大達-19.771%。其原因可能與VLW狀態(tài)方程形式的高階維里系數(shù)過度簡化有關,隨著壓力增大,描述多個氣體分子同時相互作用的高階維里系數(shù)的準確性要求提高,而VLW狀態(tài)方程中,除第二階維里系數(shù)與理論值相符合外,隨著維里系數(shù)階級的增大,其與理論值的差距也逐漸增大,故CO2的體積計算結果偏差隨壓力增大而增大。而本文所提出的對比態(tài)VHL狀態(tài)方程的第三階、第四階維里系數(shù)均與理論值吻合,高階維里系數(shù)則通過甲烷高溫高壓熱力學狀態(tài)數(shù)據(jù)優(yōu)化獲得,其具有扎實物理基礎,因此,該狀態(tài)方程能夠很好描述CO2高溫狀態(tài)下較低壓力范圍內的熱力學狀態(tài)。

表3VHL、VLW狀態(tài)方程計算CO2高溫pVT關系

Table3ThepVTrelation of CO2at high temperature calculated by the means of VHL EOS and VLW EOS

No. T/K p/GPaV/cm3·mol-1MD[12]VHLVLWerror/% VHL VLW1718.80.51164040.9118.132.275-54.6752798.80.58294040.3125.830.775-35.4253920.20.63964040.7331.431.825-21.4254963.90.68354040.3632.260.900-19.350511130.77324040.4335.221.075-11.950612300.84604040.4436.721.100-8.200712960.92084039.8536.74-0.375-8.1508720.60.77453535.9514.892.714-57.4579774.70.84463535.5320.111.514-42.54310902.00.94643535.4625.661.314-26.68611805.71.4803030.3617.631.200-41.23312999.41.7433029.9922.55-0.033-24.8331311901.9913029.8024.87-0.667-17.1001415092.2903030.0427.480.133-8.4001516752.4393030.1728.400.567-5.3331618072.6683029.8128.49-0.633-5.0331719452.8323029.7528.86-0.833-3.80018773.82.1682727.2613.790.368-49.22719887.32.3622727.0217.49-0.515-35.6042010162.4882727.1120.00-0.184-26.3622111683.6992524.7119.42-1.160-22.3202219834.9522524.7723.39-0.920-6.4402327196.0672524.7224.72-1.120-1.1202436377.0652525.0426.030.1604.1202544197.9142525.1826.720.7206.8802650608.0822525.8727.813.48011.240

Table3continued

No. T/K p/GPaV/cm3·mol-1MD[12]VHLVLWerror/% VHL VLW27700.84.4822222.662.7994.040-87.14928812.94.7212222.4511.923.076-45.27129913.44.9552222.3114.032.433-35.5833010245.2902222.1215.431.561-29.1553112187.9302020.1915.350.950-23.25032186410.522019.5617.35-2.200-13.25033262812.182019.6318.86-1.850-5.70034432115.052019.9220.70-0.4003.50035497816.362019.8821.00-0.6005.00036803.211.351818.238.881.278-50.68937122812.471818.0013.310.000-26.05638164113.561818.0015.070.000-16.27839197714.521817.9915.94-0.056-11.44440241415.641818.0216.770.111-6.83341275616.741817.9617.17-0.222-4.61142316117.561818.0417.680.222-1.77843346918.241818.0617.990.333-0.05644783.019.471616.147.190.875-55.08745118620.621615.8811.14-0.750-30.37546164521.921615.8812.91-0.750-19.31347200622.941615.9213.79-0.500-13.81348231824.111615.9014.31-0.625-10.56249277025.071616.0015.010.000-6.18750307226.321615.9415.27-0.375-4.56351338427.041615.9915.59-0.063-2.56252395627.991616.1216.160.7501.00053803.726.221515.097.010.600-53.28054119527.391514.8510.27-1.000-31.53355160828.631514.8311.77-1.133-21.53356206929.971514.8712.81-0.867-14.60057240030.921514.9113.36-0.600-10.93358271131.951514.9213.77-0.533-8.20059307532.891514.9714.19-0.200-5.40060340133.671515.0114.520.067-3.20061393834.741515.1015.010.6670.06762787.235.671414.126.170.857-55.90063121937.021413.849.51-1.143-32.07164165938.361413.8110.91-1.357-22.07165205039.581413.8311.71-1.214-16.35766234040.501413.8512.17-1.071-13.07167279341.971413.8912.74-0.786-9.00068319042.681413.9613.19-0.286-5.78669340144.181413.9013.29-0.714-5.07170789.649.451313.145.701.077-56.12371119250.771312.868.57-1.077-34.05472159452.071312.819.80-1.462-24.64673198453.321312.8210.57-1.385-18.69274231354.381312.8411.07-1.231-14.84675288656.211312.8711.74-1.000-9.69276307957.431312.8511.89-1.154-8.53877355158.001312.9312.34-0.538-5.07778395658.711312.9912.67-0.077-2.53879782.069.771212.215.111.750-57.425

Table3continued

No. T/K p/GPaV/cm3·mol-1MD[12]VHLVLWerror/% VHL VLW80120871.241211.917.86-0.750-34.52581161372.641211.848.94-1.333-25.50882199173.981211.839.61-1.417-19.91783240775.371211.8410.17-1.333-15.25084285476.901211.8610.64-1.167-11.33385323878.241211.8710.97-1.083-8.58386349879.161211.8811.16-1.000-7.00087807.0100.71111.264.962.364-54.882881209102.21110.987.12-0.182-35.264891599103.61110.908.06-0.909-26.755901971104.91110.888.66-1.091-21.273912416106.61110.879.20-1.182-16.400922866108.11110.889.63-1.091-12.500933012109.21110.879.73-1.182-11.545943326109.61110.899.98-1.000-9.236953959111.11110.9210.41-0.727-5.364

Note: error=100×(VEOS-VMD)/VMD

表41000 K時CO2pVT關系的VHL、VLW狀態(tài)方程計算值及與NIST數(shù)據(jù)庫數(shù)據(jù)的比較

Table4Comparison of the calculated values of CO2pVTrelation at 1000 K by VHL EOS and VLW EOS and the data of CO2in NIST database

No.p/GPa V/cm3·mol-1 NIST[18] VLW VHLerror/% VLW VHL10.02434.60436.9437.400.5290.64820.04228.77230.1231.500.5811.18730.06161.37161.6163.700.1431.47240.08128.35127.5130.30-0.6621.52750.10108.89107.0110.50-1.7361.43860.1296.07393.3597.31-2.8341.28970.1486.96983.5787.96-3.9081.13580.1680.15376.1980.94-4.9440.98590.1874.84270.4175.48-5.9220.847100.2070.57865.7671.08-6.8260.716110.2267.06961.9167.47-7.6920.595120.2464.12458.6864.43-8.4900.482130.2661.61355.9261.84-9.2400.374140.2859.44253.5359.60-9.9460.270150.3057.54251.4457.64-10.6040.172160.3255.86249.5855.91-11.2460.079170.3454.36447.9354.36-11.835-0.009180.3653.01746.4552.97-12.387-0.092190.3851.79845.1151.71-12.912-0.170200.4050.68743.8950.56-13.410-0.242210.4249.66942.7749.52-13.890-0.309220.4448.73241.7548.55-14.327-0.371230.4647.86640.8047.66-14.762-0.430240.4847.06239.9246.83-15.176-0.484250.5046.31239.1146.07-15.551-0.533260.5245.61138.3545.35-15.919-0.578270.5444.95437.6444.67-16.270-0.621280.5644.33636.9744.04-16.614-0.660290.5843.75436.3543.45-16.922-0.697300.6043.20435.7642.89-17.230-0.732310.6242.68335.2042.36-17.532-0.764

Table4Continued

No.p/GPaV/cm3·mol-1NIST[18]VLWVHLerror/%VLWVHL320.6442.18934.6741.85-17.822-0.794330.6641.71934.1741.38-18.095-0.822340.6841.27233.7040.92-18.347-0.848350.7040.84533.2440.49-18.619-0.872360.7240.43832.8140.08-18.863-0.896370.7440.04932.4039.68-19.099-0.920380.7639.67632.0139.30-19.322-0.941390.7839.31931.6338.94-19.555-0.963400.8038.97631.2738.59-19.771-0.983

Note: error=100×(VEOS-VNIST)/VNIST.

4 結 論

本研究提出了一種基于LJ勢能函數(shù)的對比態(tài)維里型狀態(tài)方程VHL用于描述爆轟環(huán)境下CO2的高溫高壓熱力學狀態(tài)。采用VHL狀態(tài)方程計算得CO2體積平均絕對偏差為0.971%,最大偏差為4.04%,采用VLW狀態(tài)方程計算所得平均絕對偏差20.2%,最大偏差87.149%。因此,在計算CO2高溫、中高壓熱力學狀態(tài)時,VHL狀態(tài)方程的計算準確性得到了大幅度提高。VHL狀態(tài)方程計算所得體積偏差與壓力和溫度參量無明顯的相關性; VLW狀態(tài)方程計算所得體積偏差則與溫度具有明顯的相關性,隨著溫度的升高,計算體積偏差逐漸減小。

參考文獻:

[1] Mader C L. Numerical modeling of explosives and propellants[M]. CRC press, 2007: 377-380.

[2] 吳雄, 龍新平, 何碧, 等. VLW 爆轟產物狀態(tài)方程[J]. 中國科學: B 輯, 2008, 38(12): 1129-1132.

WU Xiong,LONG Xin-ping, HE Bi, et al. VLW equation of state of detonation products[J].ScienceinChinaSeriesB:Chemistry, 2008, 38(12): 1129-1132.

[3] Fickett W, Wood W W. Tables of the Lennard-Jones and Devonshire equation of state at high temperatures and densities[J].TheJournalofChemicalPhysics, 1952, 20(10): 1624-1626.

[4] Fried L E, Howard W M. An accurate equation of state for the exponential-6 fluid applied to dense supercritical nitrogen[J].TheJournalofChemicalPhysics, 1998, 109(17): 7338-7348.

[5] 吳雄, 龍新平, 何碧, 等. VLW狀態(tài)方程的回顧與展望[J]. 高壓物理學報, 1999, 13(1): 55-58.

WU Xiong, LONG Xin-ping, HE Bi, et al. Review and look forward to the progress of VLW equation of state[J].ChineseJournalofHighPressurePhysics, 1999, 13(1): 55-58.

[6] 龍新平, 何碧, 蔣曉華, 等. 論VLW狀態(tài)方程[J]. 高壓物理學報, 2003, 17(4): 247-254.

LONG Xin-ping, HE Bi, JIANG Xiao-hua, et al. Discussions on the VLW equation of state[J].ChineseJournalofHighPressurePhysics, 2003, 17(4): 247-254.

[7] 韓勇, 龍新平, 蔣治海, 等. 用VLW狀態(tài)方程計算水的沖擊Hugoniot曲線[J]. 爆炸與沖擊, 2010, 30(1): 17-20.

HAN Yong, LONG Xin-ping, JIANG Zhi-hai, et al. A theoretical calculation for the hugoniot of water using the VLW equation of state[J].ExplosionandShockWaves, 2010, 30(1): 17-20.

[8] Barker J A, Leonard P J, Pompe A. Fifth virial coefficients[J].TheJournalofChemicalPhysics, 1966, 44(11): 4206-4211.

[9] 韓勇, 龍新平, 郭向利. 一種簡化維里型狀態(tài)方程預測高溫甲烷pVT關系[J]. 物理學報, 2014, 63(15): 150505.

HAN Yong, LONG Xin-ping, GUO Xiang-li. Prediction of methane pVT relations at high temperatures by a simplified virial equation of state[J].ActaPhysSin, 2014, 63(15): 150505.

[10] 韓勇, 龍新平, 黃毅民, 等. L-J, Exp-6 兩種形式勢能函數(shù)對計算無量綱第二維里系數(shù)的影響[J]. 含能材料, 2009, 17(5): 574-577.

HAN Yong, LONG Xin-ping, HUANG Yi-min, et al. Effect of L-J or Exp-6 potential function on calculation of reduced second Viral coefficient[J].ChineseJournalofEnergeticMaterials(HannengCailliao), 2009, 17(5): 574-577.

[11] Duan Z, M?ller N, Weare J H. Molecular dynamics simulation of PVT properties of geological fluids and a general equation of state of nonpolar and weakly polar gases up to 2000 K and 20,000 bar[J].GeochimicaetCosmochimicaActa, 1992, 56(10): 3839-3845.

[12] Belonoshko A, Saxena S K. A molecular dynamics study of the pressure-volume-temperature properties of supercritical fluids: II. CO2, CH4, CO, O2and H2[J].GeochimicaetCosmochimicaActa, 1991, 55(11): 3191-3208.

[13] Schott G L. Shock-compressed carbon dioxide: Liquid measurements and comparisons with selected models[J].HighPressureResearch, 1991, 6(3): 187-200.

[14] Nellis W J, Mitchell A C, Ree F H, et al. Equation of state of shock-compressed liquds: Carbon dioxide and air[J].JChemPhys, 1991, 95: 5268-5272.

[15] 劉福生, 陳先猛, 陳攀森, 等. 液態(tài)CO2高溫高密度狀態(tài)方程研究[J]. 高壓物理學報, 1998, 12(1): 28-33.

LIU Fu-sheng, CHEN Xian-meng, CHEN Pan-sen, et al. Equation of state of liquod CO2at highe temperatures and high densities[J].ChineseJournalofHighPressurePhysics, 1998, 12(1): 28-33.

[16] 徐士良. 常用算法程序集[M]. 北京: 清華大學出版社, 1995: 414.

[17] Ben-Amotz D, Herschbach D R. Estimation of effective diameters for molecular fluids[J].JournalofPhysicalChemistry,1990, 94(3): 1038-1047.

[18] Span R, Wagner W. A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1100 K at pressures up to 800 MPa[J].JournalofPhysicalandChemicalReferenceData, 1996, 25(6): 1509-1596.