藍(lán)文鴻， 于新國

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

符號說明：

ρ——密度，kg/m3

h——比焓，kJ/kg

T——熱力學(xué)溫度，K

w——當(dāng)?shù)芈曀伲琺/s

p——壓強(qiáng)，MPa

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

u——熱力學(xué)能，kJ/kg

v——比體積，m3/kg

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

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

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

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

1 理論推導(dǎo)及分析

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

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

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

(3)

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

2 程序設(shè)計與功能

2.1 程序設(shè)計

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

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

本文程序的主要工作原理是基于迭代式的運(yùn)算

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 程序驗證

3.1 程序計算精度驗證

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

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

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

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

根據(jù)IAPWS-IF97標(biāo)準(zhǔn)規(guī)定，滿足工業(yè)及科學(xué)計算的水和水蒸氣物性參數(shù)計算結(jié)果與標(biāo)準(zhǔn)值相對誤差應(yīng)小于0.96%。根據(jù)文獻(xiàn)[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ù)計算開源程序進(jìn)行對比來驗證本文程序的運(yùn)算速度。選用程序為哥德堡大學(xué)開發(fā)的XSteam水和水蒸氣物性參數(shù)計算程序和加州大學(xué)伯克利分校開發(fā)的IAPWS-IF97計算程序。在計算精度一致的條件下對35 000組無規(guī)律的水和水蒸氣(p,T)數(shù)據(jù)進(jìn)行計算，對3個程序進(jìn)行儲存所占空間、批量計算時占用系統(tǒng)資源和運(yùn)算速度進(jìn)行對比，結(jié)果見表5。

表5 3個程序性能指標(biāo)對比Tab.5 Performance comparison for three programs

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

從表5可以看出，本文程序在計算性能上優(yōu)于XSteam和IAPWS-IF97計算程序。本文程序儲存所占空間比XSteam計算程序小了47.02%，比IAPWS-IF97計算程序小了11.52%；本文程序的運(yùn)算時長為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)。通過適應(yīng)所測參數(shù)的需要，采用已編寫好的專門的參數(shù)轉(zhuǎn)換模塊，對程序的輸入?yún)?shù)進(jìn)行適當(dāng)調(diào)整，克服了Xsteam計算程序在非常見參數(shù)輸入下結(jié)果偏差較大的問題。因此本文程序能適應(yīng)不同場合可測參數(shù)的變化。

4 結(jié) 論

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

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

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

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