物質熱力學函數溫度系數的通用計算方法

劉 洌 梁國柱

(北京航空航天大學 宇航學院，北京 100191)

在火箭發動機的熱力特性分析中，需要使用推進劑及其燃燒產物在一定溫度下的熱力函數(摩爾定壓熱容、焓、熵)值.由于熱力學性質表大多只給出的整百K下的標準熱力函數值，當計算某一非整百K溫度的熱力函數值時，就必須利用2個或若干個相鄰溫度的標準數據進行插值，導致較大誤差產生.因此，1961年，美國的研究人員提出，考慮到物質的摩爾定壓熱容與焓、熵的相互耦合關系，可用最小二乘法擬合出計算熱力函數的多項式溫度系數，達到方便求解某一溫度范圍內，任意溫度下熱力函數的目的［1］.20世紀七八十年代，利用相似的方法，美國劉易斯研究中心通過計算機程序計算得到常用化學物質的溫度系數，并提高了計算熱力函數值的精度［2－3］.20 世紀90 年代，我國的科技工作者通過精確計算，獲得了Cu及其氧化物、氟化物的溫度系數［4］.最近20年來，美國相關研究中心進一步擴充了具有廣泛用途(包括火星大氣性質分析［5－6］)的求解物質熱力函數的溫度系數數據庫［7－10］.近年來，火箭發動機推進劑配方不斷改進，以及類似于光譜測量、原子振動測試等實驗技術的發展，不僅使火箭發動機推進劑與燃燒產物中可能存在的物質數量大幅增加，而且通過實驗測定得到的熱力函數值也更加精確，因此應該對現有物質溫度系數表進行進一步的更新和擴充.

本文根據標準熱力學數據［11－12］計算現代火箭發動機中常見的135種單質、化合物以及帶電離子在300～5 000 K溫度范圍內的溫度系數，并對計算得到的物質摩爾定壓熱容誤差較大的溫度系數進行單獨修正，為工程計算提供更準確的熱力函數溫度系數值.

1 計算方法

1.1 熱力函數的溫度系數表達

根據文獻［3］，并考慮工程中實際應用要求，可采用含7個相同溫度系數的經驗公式計算上述3種熱力函數公式，具體公式表達如下.

摩爾定壓熱容:

焓:

熵:

由熱力學微分關系形式可知，上述3種熱力函數具有關聯性，應作為整體考慮.熱力函數都具有相同的溫度系數a1～a5，焓、熵具有各自獨立的系數a6，a7;JANAF標準熱力函數表(NIST-JANAF thermochemical tables)的實驗測得值與經驗公式計算值之間具有一定的誤差.為了減少經驗公式計算值與標準值間的相對誤差，可以采用最小二乘法進行擬合計算，求出溫度系數a1～a7.

實際上，隨著溫度不斷升高，分子振動加劇，摩爾定壓熱容、焓、熵等熱力函數值發生劇烈變化，尤其要指出的是摩爾定壓熱容，對某些物質，甚至會改變摩爾定壓熱容以溫度為自變量的單調性.因此可人為將溫度分為若干區間，對于相同物質的不同的溫度區間，采用不同的溫度系數來計算熱力函數.對大多數物質而言，可分為高、低溫2個區間，對于部分固體和液體物質，需要將溫度區間更細劃分，以滿足精度要求.

1.2 溫度系數的求解方法

選取溫度區間［T0，T1］，［T1，T2］，…，［Tn－1，Tn］，共有對應的 a0，i，a1，i，…，an－1，i(i=1，2，…，7)，7n個待求溫度系數.對某種確定的物質而言，需依據不同溫度區間計算不同的溫度系數.

此外，為保證函數的精確性，要求通過相鄰溫度區間的不同溫度系數在溫度交界點計算得到的熱力函數值相等，即使分段函數兩端取得的熱力函數值相等，則必須滿足如下3(n－1)個方程:

為方便表達，用x表示式(1)～式(3)中與溫度系 數 相 乘 的 各 組 合 量，如 xCp，m，6=0，xH，2=0.5RT2，xS，1=R ln T，本文下標表示熱力函數和與之對應的溫度系數下標.則在溫度分段交接點上，利用上述熱力函數值相等的條件，可將方程變為

對式(4)～ 式(6)而言，xCp，m，i，xH，i，xS，i(i=1，2，…，7)均是交接溫度下的已知參數.

對于一個確定的j，式(4)～式(6)是關于溫度系數 aj－1，1，aj－1，2和 aj－1，3的三元一次方程組.當j=n－1時，求出其左端關于待求溫度系數的已知系數矩陣的逆矩陣B(3，3)，將其左乘于方程右端各項，獲得 an－2，k(k=1，2，3)的表達式:

由于在各溫度區間交接點上，矩陣A內均為與交接溫度相關的可計算參數.至此，如果已知式(7)右端中的4n+3個溫度系數，則由該方程即可求出左端的3(n－1)個余下的溫度系數.因此，下面將討論由求解7n個待求溫度系數變成求解4n+3個溫度系數的求解過程.

當 T ∈［Tn－1，Tn］時:

熱力函數計算值對溫度系數的偏導數為

式中xm，i僅是溫度的函數.

由式(7)，當 T ∈ ［T0，T1］∪ … ∪ ［Tn－2，Tn－1］，對應區間的溫度系數為 aj－1，i(i=1，2，…，7)時:

熱力函數計算值對溫度系數的偏導數為

由式(8)、式(9)，可得

根據偏導數關系式(11)、式(12)、式(14)和式(15)，結合式(16)、式(17)，可得到一個4n+3階的正規線性方程組，用選主元的Doolittle方法，精確地解出各個溫度區間內的4n+3個溫度系數aj，i，然后根據確定的相互關系式(7)，解出剩余的3(n－1)個未知數.

本方法能使3種熱力函數同時得到很好的擬合，又保證了在若干個不同溫度區間的交界處計算的函數值相等，誤差小，精度高.

1.3 摩爾定壓熱容的單獨修正

用1.2節中方法得到的溫度系數可以計算出對應溫度下的摩爾定壓熱容、焓和熵.但根據熱力學分析可知，作為以溫度為單自變量的函數，物質的焓、熵具有單調性，但對于某些特殊物質，摩爾定壓熱容的單調性不成立，以致使用該溫度系數計算得到的摩爾定壓熱容值誤差相對較大.其次，某些物質的摩爾定壓熱容與焓、熵數量級相差較大，也會對摩爾定壓熱容的計算精度造成影響.

用計算得到的溫度系數計算熱力函數值，與標準值進行比較，可知摩爾定壓熱容的相對誤差均比焓、熵的相對誤差大一個數量級以上，特別是在函數的單調性發生變化的時候，誤差更為明顯.為了保證數據的精確性，本文進一步對摩爾定壓熱容進行單獨修正.

為了保證擬合摩爾定壓熱容函數的精確性，可重新確定一組新的溫度系數來單獨計算物質的摩爾定壓熱容.

有最小，就可以達到修正的目的.

由式(1)可知，存在5n個待求溫度系數，根據1.2節中類似方法，通過函數變換，化簡為(4n+1)×(4n+1)的標準矩陣方程形式，從而求得待求溫度系數值.用新求解的溫度系數計算摩爾定壓熱容值，可進一步提高摩爾定壓熱容值的數據精度.

2 計算結果與分析

根據135種火箭發動機常見物質在300～5000K上的最新標準熱力函數數據，統一將溫度劃分為低溫區間(300～1 000 K)和高溫區間(1000～5000K)，采用上述計算方法，用 FORTRAN90語言編寫計算熱力函數溫度系數的程序，計算得到火箭發動機135種常見物質的溫度系數值，以及相應的摩爾定壓熱容、焓、熵的計算值，并將其與已知的熱力函數標準值進行比較，獲得上述每一種物質在各計算溫度下的熱力函數相對誤差及標準差.

現代的固體火箭發動機常采用硼和銅元素作燃燒添加劑，以B2O2與CuO氣體為例，其溫度系數計算結果如表 1、表2所示.利用式(1)～式(3)計算在上述各溫度區間上B2O2和CuO的熱力函數值，熱力函數標準值、計算值結果如圖1所示.

表1 B2 O2(氣)的溫度系數表

表2 CuO(氣)的溫度系數表

將得到的函數的計算值與標準值比較，可得如圖2所示的計算區間上各計算溫度點的相對誤差，并得到表3、表4所示的熱力函數計算值與標準值的誤差分析表.

表3 B2O2(氣)熱力函數計算值與標準值誤差分析表

表4 CuO(氣)熱力函數計算值與標準值誤差分析表

由圖2可知，對B2O2與CuO而言，摩爾定壓熱容函數值在300，1000以及5000K(2個區間端點以及交接點)上相對誤差較大，B2O2在300K時取得最大相對誤差4.13%，CuO在5000K時取得最大相對誤差3.07%，而在其他溫度下誤差均在1%以內，2種物質的平均相對誤差均在0.5%左右;而就焓、熵而言，計算值和標準值相比，在誤差較大的區間端點處相對誤差也小于1‰，計算值與標準值幾乎完全一致.

對135種物質分析計算得到的誤差結果如表5所示，可以看出利用溫度系數計算得到的計算值相對誤差均較小，完全適用于工程應用.

圖1 B2 O2(氣)、CuO(氣)熱力函數計算值與標準值比較

表5 135種物質熱力函數計算值與標準值誤差分析表

135種物質中的26種物質(氮原子、液體鉛、固體硅等)的摩爾定壓熱容計算值，從平均相對誤差而言，誤差小于5%，但最大相對誤差大于10%，有必要進行進一步修正.

圖2 熱力函數計算值與標準值相對誤差

以氮原子為例，在溫度小于1 900 K時，摩爾定壓熱容數值基本不發生變化，但隨著溫度的繼續升高，摩爾定壓熱容迅速增大.在這種條件下，若不對計算摩爾定壓熱容的溫度系數進行修正，計算中會產生很大誤差，但焓、熵的相對誤差仍然保持較低值.修正后與修正前得到的摩爾定壓熱容值相比，可以得到圖3所示結果.利用重新修正的計算摩爾定壓熱容的溫度系數計算摩爾定壓熱容，與標準摩爾定壓熱容值相比較，可以得到表6所示結果.

圖3 氮原子摩爾定壓熱容修正前后比較

表6 氮原子摩爾定壓熱容修正前后誤差比較

由表6可知，與修正前相比，整個溫度區間上摩爾定壓熱容的平均相對誤差減小了250倍，達到10－5數量級，與焓、熵精度相同，最大相對誤差也大大減小至3?，進一步提高了摩爾定壓熱容溫度系數的精確度與計算上的適應性，完全滿足工程計算需要.

當然，經驗公式計算的準確性還依賴于熱力函數標準數據的準確性.根據不同物質的特有化學性質，對單質、化合物等還可以采用3系數、9系數等各種形式的經驗公式進行熱力函數值的計算，以便進一步提高精確性.此外，根據物質不同的物理性質，劃分不同的溫度區間，也可以進一步減小計算誤差.

3 結論

為了達到方便準確計算物質熱力函數的目的，采用含7個溫度系數的多項式作為計算摩爾定壓熱容、焓、熵3種熱力函數的經驗公式，得到如下成果:

1)利用最小二乘法，計算了火箭發動機研究中常見的包括單質、化合物以及帶電離子在內的135種物質在300～5 000K溫度范圍內的溫度系數.

2)根據得到的溫度系數計算了135種物質的摩爾定壓熱容、焓、熵.與標準值相比，在任意溫度區間下，焓與熵的相對誤差均小于10－3，除極個別溫度外，摩爾定壓熱容的計算值相對誤差均小于1%.

3)對氮原子等26種摩爾定壓熱容值隨溫度變化較大的物質進行摩爾定壓熱容的溫度系數的單獨修正，修正后摩爾定壓熱容相對誤差接近10－5.

4)建立了更準確的熱力函數數據庫及溫度系數數據庫.所獲得的溫度系數不僅可運用于火箭發動機的熱力分析，還可以應用于熱能動力、化工冶金等其他學科，具有廣泛的應用價值.

References)

［1］ Zelenznik F J，Gordon S.Simultaneous least-squares approximation of a function and its first integrals with application to thermodynamic data［R］.NASA TN D-767，1961

［2］ McBride B J，Heimel S，Ehlers JG，et al.Thermodynamic properties to6000K for210 substances involving the first18 elements［R］.NASA SP-3001，1963

［3］ McBride B J，Gordon S.Computer program for calculation of complex chemical equilibrium compositions rocket performance incident and reflected shocks and Chapman-Jouguet detonations［R］.NASA SP-273，1971

［4］田德余，翁武軍，劉振源.熱力學函數溫度系數的計算［J］.航空動力學報，1991，6(2):101 －104 Tian Deyu，Weng Wujun，Liu Zhenyuan.Calculation of the thermodynamic functions’temperature coefficient［J］.Journal of Aerospace Power，1991，6(2):101 －104(in Chinese)

［5］ Captelli M，Colonna G，Giordano D.High-temperature thermodynamic properties of Mars-atmosphere components［J］.Journal of Spacecraft and Rockets，2005，42(6):980 －989

［6］ Nelson H F.Thermodynamic properties of the low Martian atmosphere［J］.Journal of Thermophysics and Heat Transfer，2005，19(4):592－594

［7］ McBride B J，Gordon S.Computer program for calculating and fitting thermodynamic functions［R］.NASA RP-1271，1992

［8］ Gordon S，McBride B J.Thermodynamic data to 20 000 K for monoatomic gases［R］.NASA TP-1999-0208523，1999

［9］ Zehe M J，Gordon S，McBride B J.CAP:a computer code for generating tabular thermodynamic functions from NASA Lewis coefficients［R］.NASA/TP-2001-210959/REV1，2001

［10］ McBride B J，Zehe M J，Gordon S.NASA glenn coefficients for calculating thermodynamic properties of individual species［R］.NASA TP-2002-211556，2002

［11］ Chase MW，Davies C A，Downey JR.JANAF thermochemical tables［M］.3 rd ed.Washington:American Chemical Society，1985:1－1856

［12］ Chase MW.NIST-JANAF thermochemical tables［M］.4 th ed.Berlin:Springer Verlag，2000:1 －1952