物質(zhì)熱力學(xué)函數(shù)溫度系數(shù)的通用計(jì)算方法

劉 洌 梁國(guó)柱

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

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

本文根據(jù)標(biāo)準(zhǔn)熱力學(xué)數(shù)據(jù)［11－12］計(jì)算現(xiàn)代火箭發(fā)動(dòng)機(jī)中常見(jiàn)的135種單質(zhì)、化合物以及帶電離子在300～5 000 K溫度范圍內(nèi)的溫度系數(shù)，并對(duì)計(jì)算得到的物質(zhì)摩爾定壓熱容誤差較大的溫度系數(shù)進(jìn)行單獨(dú)修正，為工程計(jì)算提供更準(zhǔn)確的熱力函數(shù)溫度系數(shù)值.

1 計(jì)算方法

1.1 熱力函數(shù)的溫度系數(shù)表達(dá)

根據(jù)文獻(xiàn)［3］，并考慮工程中實(shí)際應(yīng)用要求，可采用含7個(gè)相同溫度系數(shù)的經(jīng)驗(yàn)公式計(jì)算上述3種熱力函數(shù)公式，具體公式表達(dá)如下.

摩爾定壓熱容:

焓:

熵:

由熱力學(xué)微分關(guān)系形式可知，上述3種熱力函數(shù)具有關(guān)聯(lián)性，應(yīng)作為整體考慮.熱力函數(shù)都具有相同的溫度系數(shù)a1～a5，焓、熵具有各自獨(dú)立的系數(shù)a6，a7;JANAF標(biāo)準(zhǔn)熱力函數(shù)表(NIST-JANAF thermochemical tables)的實(shí)驗(yàn)測(cè)得值與經(jīng)驗(yàn)公式計(jì)算值之間具有一定的誤差.為了減少經(jīng)驗(yàn)公式計(jì)算值與標(biāo)準(zhǔn)值間的相對(duì)誤差，可以采用最小二乘法進(jìn)行擬合計(jì)算，求出溫度系數(shù)a1～a7.

實(shí)際上，隨著溫度不斷升高，分子振動(dòng)加劇，摩爾定壓熱容、焓、熵等熱力函數(shù)值發(fā)生劇烈變化，尤其要指出的是摩爾定壓熱容，對(duì)某些物質(zhì)，甚至?xí)淖兡柖▔簾崛菀詼囟葹樽宰兞康膯握{(diào)性.因此可人為將溫度分為若干區(qū)間，對(duì)于相同物質(zhì)的不同的溫度區(qū)間，采用不同的溫度系數(shù)來(lái)計(jì)算熱力函數(shù).對(duì)大多數(shù)物質(zhì)而言，可分為高、低溫2個(gè)區(qū)間，對(duì)于部分固體和液體物質(zhì)，需要將溫度區(qū)間更細(xì)劃分，以滿足精度要求.

1.2 溫度系數(shù)的求解方法

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

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

為方便表達(dá)，用x表示式(1)～式(3)中與溫度系 數(shù) 相 乘 的 各 組 合 量，如 xCp，m，6=0，xH，2=0.5RT2，xS，1=R ln T，本文下標(biāo)表示熱力函數(shù)和與之對(duì)應(yīng)的溫度系數(shù)下標(biāo).則在溫度分段交接點(diǎn)上，利用上述熱力函數(shù)值相等的條件，可將方程變?yōu)?/p>

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

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

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

當(dāng) T ∈［Tn－1，Tn］時(shí):

熱力函數(shù)計(jì)算值對(duì)溫度系數(shù)的偏導(dǎo)數(shù)為

式中xm，i僅是溫度的函數(shù).

由式(7)，當(dāng) T ∈ ［T0，T1］∪ … ∪ ［Tn－2，Tn－1］，對(duì)應(yīng)區(qū)間的溫度系數(shù)為 aj－1，i(i=1，2，…，7)時(shí):

熱力函數(shù)計(jì)算值對(duì)溫度系數(shù)的偏導(dǎo)數(shù)為

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

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

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

1.3 摩爾定壓熱容的單獨(dú)修正

用1.2節(jié)中方法得到的溫度系數(shù)可以計(jì)算出對(duì)應(yīng)溫度下的摩爾定壓熱容、焓和熵.但根據(jù)熱力學(xué)分析可知，作為以溫度為單自變量的函數(shù)，物質(zhì)的焓、熵具有單調(diào)性，但對(duì)于某些特殊物質(zhì)，摩爾定壓熱容的單調(diào)性不成立，以致使用該溫度系數(shù)計(jì)算得到的摩爾定壓熱容值誤差相對(duì)較大.其次，某些物質(zhì)的摩爾定壓熱容與焓、熵?cái)?shù)量級(jí)相差較大，也會(huì)對(duì)摩爾定壓熱容的計(jì)算精度造成影響.

用計(jì)算得到的溫度系數(shù)計(jì)算熱力函數(shù)值，與標(biāo)準(zhǔn)值進(jìn)行比較，可知摩爾定壓熱容的相對(duì)誤差均比焓、熵的相對(duì)誤差大一個(gè)數(shù)量級(jí)以上，特別是在函數(shù)的單調(diào)性發(fā)生變化的時(shí)候，誤差更為明顯.為了保證數(shù)據(jù)的精確性，本文進(jìn)一步對(duì)摩爾定壓熱容進(jìn)行單獨(dú)修正.

為了保證擬合摩爾定壓熱容函數(shù)的精確性，可重新確定一組新的溫度系數(shù)來(lái)單獨(dú)計(jì)算物質(zhì)的摩爾定壓熱容.

有最小，就可以達(dá)到修正的目的.

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

2 計(jì)算結(jié)果與分析

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

現(xiàn)代的固體火箭發(fā)動(dòng)機(jī)常采用硼和銅元素作燃燒添加劑，以B2O2與CuO氣體為例，其溫度系數(shù)計(jì)算結(jié)果如表 1、表2所示.利用式(1)～式(3)計(jì)算在上述各溫度區(qū)間上B2O2和CuO的熱力函數(shù)值，熱力函數(shù)標(biāo)準(zhǔn)值、計(jì)算值結(jié)果如圖1所示.

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

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

將得到的函數(shù)的計(jì)算值與標(biāo)準(zhǔn)值比較，可得如圖2所示的計(jì)算區(qū)間上各計(jì)算溫度點(diǎn)的相對(duì)誤差，并得到表3、表4所示的熱力函數(shù)計(jì)算值與標(biāo)準(zhǔn)值的誤差分析表.

表3 B2O2(氣)熱力函數(shù)計(jì)算值與標(biāo)準(zhǔn)值誤差分析表

表4 CuO(氣)熱力函數(shù)計(jì)算值與標(biāo)準(zhǔn)值誤差分析表

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

對(duì)135種物質(zhì)分析計(jì)算得到的誤差結(jié)果如表5所示，可以看出利用溫度系數(shù)計(jì)算得到的計(jì)算值相對(duì)誤差均較小，完全適用于工程應(yīng)用.

圖1 B2 O2(氣)、CuO(氣)熱力函數(shù)計(jì)算值與標(biāo)準(zhǔn)值比較

表5 135種物質(zhì)熱力函數(shù)計(jì)算值與標(biāo)準(zhǔn)值誤差分析表

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

圖2 熱力函數(shù)計(jì)算值與標(biāo)準(zhǔn)值相對(duì)誤差

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

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

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

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

當(dāng)然，經(jīng)驗(yàn)公式計(jì)算的準(zhǔn)確性還依賴于熱力函數(shù)標(biāo)準(zhǔn)數(shù)據(jù)的準(zhǔn)確性.根據(jù)不同物質(zhì)的特有化學(xué)性質(zhì)，對(duì)單質(zhì)、化合物等還可以采用3系數(shù)、9系數(shù)等各種形式的經(jīng)驗(yàn)公式進(jìn)行熱力函數(shù)值的計(jì)算，以便進(jìn)一步提高精確性.此外，根據(jù)物質(zhì)不同的物理性質(zhì)，劃分不同的溫度區(qū)間，也可以進(jìn)一步減小計(jì)算誤差.

3 結(jié)論

為了達(dá)到方便準(zhǔn)確計(jì)算物質(zhì)熱力函數(shù)的目的，采用含7個(gè)溫度系數(shù)的多項(xiàng)式作為計(jì)算摩爾定壓熱容、焓、熵3種熱力函數(shù)的經(jīng)驗(yàn)公式，得到如下成果:

1)利用最小二乘法，計(jì)算了火箭發(fā)動(dòng)機(jī)研究中常見(jiàn)的包括單質(zhì)、化合物以及帶電離子在內(nèi)的135種物質(zhì)在300～5 000K溫度范圍內(nèi)的溫度系數(shù).

2)根據(jù)得到的溫度系數(shù)計(jì)算了135種物質(zhì)的摩爾定壓熱容、焓、熵.與標(biāo)準(zhǔn)值相比，在任意溫度區(qū)間下，焓與熵的相對(duì)誤差均小于10－3，除極個(gè)別溫度外，摩爾定壓熱容的計(jì)算值相對(duì)誤差均小于1%.

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

4)建立了更準(zhǔn)確的熱力函數(shù)數(shù)據(jù)庫(kù)及溫度系數(shù)數(shù)據(jù)庫(kù).所獲得的溫度系數(shù)不僅可運(yùn)用于火箭發(fā)動(dòng)機(jī)的熱力分析，還可以應(yīng)用于熱能動(dòng)力、化工冶金等其他學(xué)科，具有廣泛的應(yīng)用價(jià)值.

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］田德余，翁武軍，劉振源.熱力學(xué)函數(shù)溫度系數(shù)的計(jì)算［J］.航空動(dòng)力學(xué)報(bào)，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