水和水蒸氣物性參數計算程序開發

藍文鴻， 于新國

(華北電力大學 核科學與工程學院，北京 102206)

符號說明：

ρ——密度，kg/m3

h——比焓，kJ/kg

T——熱力學溫度，K

w——當地聲速，m/s

p——壓強，MPa

s——比熵，kJ/(kg·K)

u——熱力學能，kJ/kg

v——比體積，m3/kg

cp——比定壓熱容，kJ/(kg·K)

cv——比定容熱容，kJ/(kg·K)

精確快速計算出水和水蒸氣的物性參數對核電站的設計和熱工水力分析程序開發至關重要。許多高校和研究所均有自己獨立開發的水和水蒸氣物性參數計算程序[1]，這些程序大多基于C++和Fortran語言編寫[2-3]，能對少量數據進行精確計算，精度要求滿足目前工業應用需要。但現有的IAPWS-IF97標準對水和水蒸氣物性參數的計算主要側重于優化迭代速度，降低迭代次數[4]，較少涉及到因程序截斷誤差引起的計算結果不準確的分析。同時現有程序也沒有對計算結果進行校核及進一步的檢驗，導致由于精度要求的提高，其截斷誤差隨著迭代次數的增加而積累，從而引起結果不準確的問題。此外，目前的程序大多采用手動少量數據輸入進行運算的模式，不具備對大規模的數據進行批量導入、計算和導出的功能。

過去基于工程語言編寫的程序存在計算精度和運算速度不足的問題，而Matlab語言簡單易懂，能夠以圖形化界面直觀展示大規模計算結果的數據變化規律，同時其對截斷誤差的處理優于C++和Fortran。在程序結構上，筆者設計了更合理的計算步驟和程序邏輯，基于導出式分析解與迭代式數值解相互校核檢驗結果的思路提出反饋校核函數，保證程序在不同精度要求下仍能穩定快速輸出計算結果。在數據輸入、輸出上設計了相關接口，能快速輸出多種數據格式，為大規模計算提供可能。水和水蒸氣物性參數計算程序基于IAPWS-IF97標準[5]進行編寫，區域劃分采用IAPWS-IF97標準推薦劃分方式，其計算參數覆蓋范圍為0≤p≤100 MPa、273.15 K≤T≤2 273.15 K。

1 理論推導及分析

根據IAPWS-IF97標準，結合文獻[6]給定的單位質量工質的亥姆霍茲自由能f=u-Ts和單位質量工質的吉布斯自由能g=h-Ts，確定水和水蒸氣狀態，依據水和水蒸氣的狀態將參數待計算范圍劃分為5個區域[5](見圖1)，分別對應：常規水區、過熱蒸汽區、臨界區、飽和水線和高溫低壓區。在臨界區，依據臨界水的狀態劃分了26個二級區域[7-8]進行計算。采用吉布斯自由能函數g(p,T)和亥姆霍茲自由能函數f(ρ,T)為基本方程，具體公式如下:

f(ρ,T)=u(ρ,T)-T×s(ρ,T)

(1)

g(p,T)=h(p,T)-T×s(p,T)

(2)

通過g(p,T)和f(ρ,T)的導出式即可計算基本的水和水蒸氣物性參數。

圖1 區域和對應方程Fig.1 Division of regions and corresponding equations of IAPWS-IF97

對于簡單可壓縮的純物質系統，任意一個狀態參數都可由另外2個相互獨立的狀態參數的特征函數表示。通過對特征函數應用全微分得到麥克斯韋關系式：

(3)

根據麥克斯韋關系式和拉普拉斯方程推導得出各物性參數，結果見表1和表2。

2 程序設計與功能

2.1 程序設計

采用偏導導出式計算時，需要逐次迭代解高階偏微分方程，雖然計算結果精確，但計算時間過長，不利于程序的快速響應。因此，本文程序沒有采用偏導導出式的解法，而是采用自擬合的迭代式的數值解法。對IAPWS-IF97標準中的擬合迭代式進行修正，得出本程序各區域的最終迭代式。

表1 g(p,T)推導得出的物性參數Tab.1 Physical properties obtained by partial derivative of g(p,T)

本文程序的主要工作原理是基于迭代式的運算

表2 f(ρ,T)推導得出的物性參數Tab.2 Physical properties obtained by partial derivative of f(ρ,T)

求得水和水蒸氣物性參數。主程序的流程圖見圖2。

圖2 程序流程圖Fig.2 Program flow chart

程序輸入的獨立參數為(p,T)，當程序讀入輸入參數(p,T)時，會先利用區域判別算法對參數所在區域進行判別,再將待求參數輸入到對應子區域進行運算。當程序在Matlab軟件下運行時，激活反饋校核算法，對迭代式計算結果與偏導導出式計算結果進行比較后輸出計算結果。當需要快速響應或不在Matlab軟件下運行時，將不會使用反饋校核算法。

反饋校核算法流程圖如圖3所示，主要用于對迭代式計算結果進行校核檢驗。反饋校核算法首先根據區域判別算法給定的區域選擇相應的吉布斯自由能函數g(p,T)或亥姆霍茲自由能函數f(ρ,T)。進而對g(p,T)或f(ρ,T)進行偏導導出式運算，求得各參數的近似分析解。通過反饋校核算法對偏導導出式計算結果與迭代式計算結果進行比較，當兩計算結果差值小于1%時，判定結果可信，輸出主程序結果；當其計算結果差值大于等于1%時，反饋校核算法將要求主程序輸入參數精度提高三位，再次運行主程序和反饋校核算法，直到計算結果差值小于1%為止。

圖3 反饋校核算法流程圖Fig.3 Feedback calculation chart

由于IAPWS-IF97在區域3的計算十分繁瑣，區域劃分多，迭代次數高，程序中并沒有沿用其他常用的遍歷區域判別算法和多次迭代的參數計算方法。筆者采用數學建模中常用的求最優解的模擬退火算法[9]中的區域判別算法函數，其算法流程圖見圖4。

圖4 區域判別算法流程圖Fig.4 Flow chart of regional discriminating algorithm

對區域判別算法函數在區域3的應用進行特殊處理。先通過轉換函數將輸入參數(ρ,T)轉化為(p,T)，將(p,T)視作坐標點，即p為橫坐標，T為縱坐標，再計算坐標點(p,T)與參考點(0，273.15)的距離L1，比較參考點到各區域中心點的距離L2與L1的大小，當L2與L1接近時(其差值小于區域長度的一半)，則認為待求參數落在該區域內。將輸入參數(p,T)與匹配區域的邊界進行比較、校核，滿足則確認待求參數落在該區域內，若不滿足則將區域移至相鄰區域再進行判別，直至找到滿足區域或判斷其超出計算邊界為止。

此外，區域判別算法在待計算參數的全局區域也能使用，經過多組數據的測試，其在全局區域判定中無需遍歷5個區域即可準確判斷出待求參數所在區域，避免了IAPWS-IF97標準中提供的全區遍歷算法須全區遍歷的缺點。

為了適應不同輸入參數的計算需求，本文程序還附帶了不同獨立參數之間的轉換模塊，各模塊的導出式也可由麥克斯韋關系式得出，因此能夠靈活適配各種輸入情況。各模塊如表3所示。

表3 輸入參數的轉換模塊Tab.3 Converting module for input variables

2.2 程序主要實現功能

編程語言采用科學計算語言Matlab。根據核電站的實際情況，程序的計算精度設定為64位，其顯示精度設定為10位。該程序主要實現以下功能：

(1)計算工質的飽和狀態參數，即根據工質的飽和溫度求取飽和壓力，或根據飽和壓力求取飽和溫度。此功能采用IAPWS-IF97標準提供的飽和壓力psat(T)和飽和溫度Tsat(p)函數實現。

(2)根據工質的溫度和壓力求取工質的物性參數。基于IAPWS-IF97標準提供的g(p,T)函數與自擬合的迭代式計算出區域1、區域2和區域5工質的物性參數。對于區域3的求解，采用國際水和水蒸氣物性協會(IAPWS)發布的IF97-S05論文[7]提供的補充方程v(p,T)先行求解工質密度ρ，再代入f(ρ,T)進行計算。

(3)程序根據實際要求，已經預留其與其他相關程序連接的接口。目前，本文程序已經能夠實現基于Aduino單片機平臺進行水和水蒸氣物性參數的計算，即本文程序能夠作為其他更大規模的分析計算軟件或者嵌入式系統，提供水和水蒸氣物性參數的計算數據。

3 程序驗證

3.1 程序計算精度驗證

選取IAPWS-IF97標準提供的推薦檢驗數值進行計算，本文程序計算結果與IAPWS-IF97標準對比數據見表4。由表4可知，通過對推薦檢驗數值的計算可以發現，所有計算結果的誤差均低于0.1%，這說明本文程序達到了IAPWS-IF97標準的精度要求。

另外，通過對比本文程序求取的300 K時不同壓力下水和水蒸氣比體積的結果與2008年版國際水和水蒸氣表[10]中相應條件下的比體積檢驗數據，來驗證所開發的程序。以上述數據為例進行誤差分析，結果如圖5所示，計算結果與標準值之間相對誤差低于文獻[11]給出的工程標準10-6%。

表4 計算結果與標準值的對比Tab.4 Comparison of calculation results with standard values

圖5 300 K時比體積標準值與計算結果的相對誤差Fig.5 Relative error of calculated specific volume over standard value at 300 K

根據IAPWS-IF97標準規定，滿足工業及科學計算的水和水蒸氣物性參數計算結果與標準值相對誤差應小于0.96%。根據文獻[11]，中廣核大亞灣嶺澳核電站要求水和水蒸氣物性參數計算相對誤差不高于0.96%，一回路關鍵物性參數(溫度T、壓力p、比體積v和比定壓熱容cp)計算相對誤差不高于0.81%。與2008年版國際水和水蒸氣表對比，本文程序計算結果相對誤差小于0.96%，關鍵物性參數計算結果的相對誤差小于0.7%。滿足一般工業及壓水堆核電廠的精度要求。

3.2 程序性能驗證

所開發的水和水蒸氣物性參數計算程序采用了反饋校核算法，在區域判別上借鑒了模擬退火算法，具有優秀的計算精度以及速度、高度的可移植性和拓展性。

通過與一些水和水蒸氣物性參數計算開源程序進行對比來驗證本文程序的運算速度。選用程序為哥德堡大學開發的XSteam水和水蒸氣物性參數計算程序和加州大學伯克利分校開發的IAPWS-IF97計算程序。在計算精度一致的條件下對35 000組無規律的水和水蒸氣(p,T)數據進行計算，對3個程序進行儲存所占空間、批量計算時占用系統資源和運算速度進行對比，結果見表5。

表5 3個程序性能指標對比Tab.5 Performance comparison for three programs

本次對比測試采用控制變量法，結合Matlab軟件自有的在線監測系統記錄其運行狀態從而得出相關對比結果。

從表5可以看出，本文程序在計算性能上優于XSteam和IAPWS-IF97計算程序。本文程序儲存所占空間比XSteam計算程序小了47.02%，比IAPWS-IF97計算程序小了11.52%；本文程序的運算時長為19.57 s，比XSteam計算程序快43.01%，比IAPWS-IF97計算程序快26.29%；對于計算時占用系統資源，本文程序峰值為18.65%，低于XSteam計算程序的23.41%，更低于IAPWS-IF97計算程序的42.81%。

在程序結構上，本文程序采用了積木式的模塊化結構。通過適應所測參數的需要，采用已編寫好的專門的參數轉換模塊，對程序的輸入參數進行適當調整，克服了Xsteam計算程序在非常見參數輸入下結果偏差較大的問題。因此本文程序能適應不同場合可測參數的變化。

4 結 論

(1)基于科學計算語言Matlab進行數值計算時精度要高于C++和Fortran，同時能夠更快、更精確地求取水和水蒸氣物性參數，程序求取結果對比2008年版國際水和水蒸氣表，其相對誤差最大不超過0.96%，符合IAPWS-IF97標準要求。求取精度為64位精度，達到輕水堆核電站熱工計算的應用要求。

(2)目前基于IAPWS-IF97標準的物性參數計算軟件，其大量機時被消耗在參數的區域判定上。本文程序借鑒模擬退火算法提出了一種區域判別算法，能實現更快、更高效的輸入響應，為實現大量數據計算和在線監測的即時反饋提供可能。

(3)在核心計算模塊上采用了盡量簡單和通用的函數，故程序易于移植，能夠在不同平臺和設備上運行。同時其所占空間很小，能夠移植至嵌入式系統進行運算。

(4)本文程序首次引入反饋校驗算法，通過比對分析解與迭代式數值解的方式，相互校核計算結果，能夠有效保證計算精度，進一步確保計算結果的可靠性。