水和水蒸氣物性參數(shù)計(jì)算程序開發(fā)

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

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

符號(hào)說(shuō)明：

ρ——密度，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)

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

過(guò)去基于工程語(yǔ)言編寫的程序存在計(jì)算精度和運(yùn)算速度不足的問(wèn)題，而Matlab語(yǔ)言簡(jiǎn)單易懂，能夠以圖形化界面直觀展示大規(guī)模計(jì)算結(jié)果的數(shù)據(jù)變化規(guī)律，同時(shí)其對(duì)截?cái)嗾`差的處理優(yōu)于C++和Fortran。在程序結(jié)構(gòu)上，筆者設(shè)計(jì)了更合理的計(jì)算步驟和程序邏輯，基于導(dǎo)出式分析解與迭代式數(shù)值解相互校核檢驗(yàn)結(jié)果的思路提出反饋校核函數(shù)，保證程序在不同精度要求下仍能穩(wěn)定快速輸出計(jì)算結(jié)果。在數(shù)據(jù)輸入、輸出上設(shè)計(jì)了相關(guān)接口，能快速輸出多種數(shù)據(jù)格式，為大規(guī)模計(jì)算提供可能。水和水蒸氣物性參數(shù)計(jì)算程序基于IAPWS-IF97標(biāo)準(zhǔn)[5]進(jìn)行編寫，區(qū)域劃分采用IAPWS-IF97標(biāo)準(zhǔn)推薦劃分方式，其計(jì)算參數(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ù)待計(jì)算范圍劃分為5個(gè)區(qū)域[5](見圖1)，分別對(duì)應(yīng)：常規(guī)水區(qū)、過(guò)熱蒸汽區(qū)、臨界區(qū)、飽和水線和高溫低壓區(qū)。在臨界區(qū)，依據(jù)臨界水的狀態(tài)劃分了26個(gè)二級(jí)區(qū)域[7-8]進(jìn)行計(jì)算。采用吉布斯自由能函數(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)

通過(guò)g(p,T)和f(ρ,T)的導(dǎo)出式即可計(jì)算基本的水和水蒸氣物性參數(shù)。

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

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

(3)

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

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

2.1 程序設(shè)計(jì)

采用偏導(dǎo)導(dǎo)出式計(jì)算時(shí)，需要逐次迭代解高階偏微分方程，雖然計(jì)算結(jié)果精確，但計(jì)算時(shí)間過(guò)長(zhǎng)，不利于程序的快速響應(yīng)。因此，本文程序沒(méi)有采用偏導(dǎo)導(dǎo)出式的解法，而是采用自擬合的迭代式的數(shù)值解法。對(duì)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)時(shí)，會(huì)先利用區(qū)域判別算法對(duì)參數(shù)所在區(qū)域進(jìn)行判別,再將待求參數(shù)輸入到對(duì)應(yīng)子區(qū)域進(jìn)行運(yùn)算。當(dāng)程序在Matlab軟件下運(yùn)行時(shí)，激活反饋校核算法，對(duì)迭代式計(jì)算結(jié)果與偏導(dǎo)導(dǎo)出式計(jì)算結(jié)果進(jìn)行比較后輸出計(jì)算結(jié)果。當(dāng)需要快速響應(yīng)或不在Matlab軟件下運(yùn)行時(shí)，將不會(huì)使用反饋校核算法。

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

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

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

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

對(duì)區(qū)域判別算法函數(shù)在區(qū)域3的應(yīng)用進(jìn)行特殊處理。先通過(guò)轉(zhuǎn)換函數(shù)將輸入?yún)?shù)(ρ,T)轉(zhuǎn)化為(p,T)，將(p,T)視作坐標(biāo)點(diǎn)，即p為橫坐標(biāo)，T為縱坐標(biāo)，再計(jì)算坐標(biāo)點(diǎn)(p,T)與參考點(diǎn)(0，273.15)的距離L1，比較參考點(diǎn)到各區(qū)域中心點(diǎn)的距離L2與L1的大小，當(dāng)L2與L1接近時(shí)(其差值小于區(qū)域長(zhǎng)度的一半)，則認(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ū)域或判斷其超出計(jì)算邊界為止。

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

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

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

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

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

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

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

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

3 程序驗(yàn)證

3.1 程序計(jì)算精度驗(yàn)證

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

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

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

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

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

3.2 程序性能驗(yàn)證

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

通過(guò)與一些水和水蒸氣物性參數(shù)計(jì)算開源程序進(jìn)行對(duì)比來(lái)驗(yàn)證本文程序的運(yùn)算速度。選用程序?yàn)楦绲卤ご髮W(xué)開發(fā)的XSteam水和水蒸氣物性參數(shù)計(jì)算程序和加州大學(xué)伯克利分校開發(fā)的IAPWS-IF97計(jì)算程序。在計(jì)算精度一致的條件下對(duì)35 000組無(wú)規(guī)律的水和水蒸氣(p,T)數(shù)據(jù)進(jìn)行計(jì)算，對(duì)3個(gè)程序進(jìn)行儲(chǔ)存所占空間、批量計(jì)算時(shí)占用系統(tǒng)資源和運(yùn)算速度進(jìn)行對(duì)比，結(jié)果見表5。

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

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

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

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

4 結(jié) 論

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

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

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

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