藍文鴻， 于新國

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

符號說明：

ρ——密度，kg/m3

h——比焓，kJ/kg

T——熱力學溫度，K

w——當?shù)芈曀?，m/s

p——壓強，MPa

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

u——熱力學能，kJ/kg

v——比體積，m3/kg

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

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

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

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

1 理論推導及分析

根據(jù)IAPWS-IF97標準，結(jié)合文獻[6]給定的單位質(zhì)量工質(zhì)的亥姆霍茲自由能f=u-Ts和單位質(zhì)量工質(zhì)的吉布斯自由能g=h-Ts，確定水和水蒸氣狀態(tài)，依據(jù)水和水蒸氣的狀態(tài)將參數(shù)待計算范圍劃分為5個區(qū)域[5](見圖1)，分別對應：常規(guī)水區(qū)、過熱蒸汽區(qū)、臨界區(qū)、飽和水線和高溫低壓區(qū)。在臨界區(qū)，依據(jù)臨界水的狀態(tài)劃分了26個二級區(qū)域[7-8]進行計算。采用吉布斯自由能函數(shù)g(p,T)和亥姆霍茲自由能函數(shù)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)的導出式即可計算基本的水和水蒸氣物性參數(shù)。

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

對于簡單可壓縮的純物質(zhì)系統(tǒng)，任意一個狀態(tài)參數(shù)都可由另外2個相互獨立的狀態(tài)參數(shù)的特征函數(shù)表示。通過對特征函數(shù)應用全微分得到麥克斯韋關(guān)系式：

(3)

根據(jù)麥克斯韋關(guān)系式和拉普拉斯方程推導得出各物性參數(shù)，結(jié)果見表1和表2。

2 程序設計與功能

2.1 程序設計

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表3 輸入?yún)?shù)的轉(zhuǎn)換模塊Tab.3 Converting module for input variables

2.2 程序主要實現(xiàn)功能

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

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

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

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

3 程序驗證

3.1 程序計算精度驗證

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

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

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

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

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

3.2 程序性能驗證

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

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

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

本次對比測試采用控制變量法，結(jié)合Matlab軟件自有的在線監(jiān)測系統(tǒng)記錄其運行狀態(tài)從而得出相關(guān)對比結(jié)果。

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

在程序結(jié)構(gòu)上，本文程序采用了積木式的模塊化結(jié)構(gòu)。通過適應所測參數(shù)的需要，采用已編寫好的專門的參數(shù)轉(zhuǎn)換模塊，對程序的輸入?yún)?shù)進行適當調(diào)整，克服了Xsteam計算程序在非常見參數(shù)輸入下結(jié)果偏差較大的問題。因此本文程序能適應不同場合可測參數(shù)的變化。

4 結(jié) 論

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

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

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

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