張銀星， 高璞珍， 何曉強， 陳沖， 林宇琦， 劉怡雯

(1.哈爾濱工程大學 核安全與仿真技術國防重點學科實驗室，黑龍江 哈爾濱 150001； 2.海軍工程大學 核科學技術學院，湖北 武漢 430033； 3.中國艦船研究設計中心，湖北 武漢 430064； 4.武漢第二船舶設計研究所，湖北 武漢 430064)

核能的發(fā)展與應用對于我國綜合國力的提升有重要意義。民用核能主要應用于核反應堆中，對于大多數壓水堆，堆芯為棒束通道結構。因此研究棒束通道內的熱工水力特性對于徹底了解核反應堆的運行機理至關重要。若反應堆在運行過程中發(fā)生了事故，非能動安全系統(tǒng)就會發(fā)揮作用以確保反應堆的安全，因此棒束通道內的自然循環(huán)流動特性是值得研究的。由于棒束通道結構的特殊性，有許多學者通過實驗、數值模擬進行研究[1-5]。實驗研究可以得到更為準確、更令人信服的實驗結果，數值模擬研究可以節(jié)省大量的人力、物力、財力。但對于想要較為精準地通過數值模擬來研究自然循環(huán)條件下棒束通道的熱工水力特性，就需要將包括棒束通道的整個自然循環(huán)回路都進行模擬，這對于三維CFD數值模擬來說需要的時間更長，對計算機的要求更高。為解決這一困難，許多學者對三維CFD數值模擬軟件與一維程序耦合進行了研究，李偉等[6]證明了三維CFD數值模擬軟件FLUENT與熱工水力分析系統(tǒng)程序RELAP5耦合研究可以很好地分析簡單的單相流或復雜的兩相流問題；Palazzi等[7-8]通過模擬管道中的摩擦阻力系數證明了STAR-CCM+與RELAP5-3D耦合的可行性；Gilman[9]通過編寫用戶自定義程序建立了新沸騰模型，該模型能夠準確預測CFD模擬中過冷沸騰流體的溫度和熱流密度。王寧波等[10]對5×5棒束通道內的壓降特性進行了數值研究，但對于單一的CFD數值模擬，僅能模擬強迫循環(huán)條件下的流動特性，無法模擬自然循環(huán)條件下的熱工水力特性。孔松等[11]用一維分析程序RELAP5對小型核動力裝置進行了自然循環(huán)運行特性分析，但由于一維程序的局限性，很難對結構復雜的部件(例如反應堆堆芯)進行分析。

為探究反應堆堆芯內的固有安全性，本文研究自然循環(huán)條件下棒束通道內的流動特性，將研究棒束通道的三維CFD數值模擬與除去棒束通道外的自然循環(huán)回路其他部分的一維程序模擬耦合計算，以期用經濟的方法計算得到更為準確的數值模擬結果。通過這種方法既可以將棒束通道通過三維CFD數值計算較為精細的模擬出來，又可以考慮自然循環(huán)回路對棒束通道的影響。自然循環(huán)回路管道多為普通圓管，采用一維程序進行模擬也可以得到很好的計算結果。本文采用三維CFD數值模擬軟件STAR-CCM+與一維用戶程序(User Code)耦合以實現(xiàn)自然循環(huán)條件下棒束通道的數值模擬。

1 三維與一維耦合方法

為確保三維與一維程序耦合數值模擬結果的準確性，首先要進行實驗研究，將得到的自然循環(huán)實驗數據用來驗證數值模擬結果準確與否。實驗裝置為包含3×3棒束通道的自然循環(huán)回路，如圖1所示。本文所述三維與一維程序耦合模擬皆基于該試驗臺。具體實驗方法見文獻[2]。

圖1 自然循環(huán)實驗裝置Fig.1 Natural circulation experimental facility

為模擬自然循環(huán)條件下棒束通道的流動特性，本文采用基于有限體積法的CFD軟件STAR-CCM+ 13.04對實驗段3×3棒束通道進行三維數值模擬，并利用編程語言C++自行編寫除去棒束通道外的自然循環(huán)回路其他管道，再將2個部分通過STAR-CCM+的用戶程序自定義接口交互，以實現(xiàn)自然循環(huán)回路的完整數值模擬。

1.1 三維CFD數值模擬

棒束通道實驗段結構如圖2所示。棒束通道實驗段長度為800 mm，利用STAR-CCM+內的3D-CAD模型可直接進行建模。

1.1.1 網格獨立性分析

本文采用多面體網格實現(xiàn)棒束通道主體結構的網格劃分，采用棱柱層網格模型應用在所有的壁面邊界。圖3給出棒束通道內生成的三維網格。

為了驗證數值計算的準確性，首先需進行網格無關性驗證。在網格獨立的條件下對數值模型進行驗證。以入口速度v=2 m/s的工況為例，對棒束通道進行了網格考核。以棒束通道出口溫度To、棒束通道出口速度vo及棒束通道進出口壓降ΔP作為監(jiān)測量。選取網格數分別為868 147、1 192 331、1 213 874和1 671 941的4套網格，當網格數從868 147變化到1 671 941時，棒束通道的出口溫度To、出口速度vo和進出口壓降ΔP隨網格數量的變化趨勢如表1所示。從表1中可以看出，模型2中To的改變量為0.003 3%，vo的改變量為0.15%，ΔP的改變量為1.03%，故認為網格數在1 192 331已達到獨立。本文計算中選取網格數為1 192 331。

圖2 棒束通道結構Fig.2 Rod bundle channel structure

圖3 棒束通道網格示意Fig.3 Rod bundle channel grid diagram

表1 網格無關性驗證Table 1 Grid independence verification

1.1.2 模型驗證

為保證數值模擬的準確性，本文對豎直條件下的棒束通道進行模型驗證。在圖1所示的實驗裝置中進行絕對壓力為0.4 MPa，棒束通道加熱功率為18 kW/m2的實驗，將得到的實驗結果與相同條件下不同湍流模型的計算結果進行比較。為找到最適合棒束通道的湍流模型，本文嘗試了6種模型，分別為標準(Wilcox)k-ω模型、標準k-ε模型、可實現(xiàn)的k-ε模型、可實現(xiàn)的k-ε兩層模型、雷諾應力湍流模型與SST (Menter)k-ω模型。數值模擬過程中將棒束通道入口設為速度入口邊界，出口設為壓力出口邊界，棒束壁面設為定熱流量，外壁面設為絕熱，由此得到棒束通道進出口壓降ΔP隨雷諾數Re的變化曲線如圖4所示。從圖4中可以看出SST (Menter)k-ω模型可以很好地模擬棒束通道內的壓降特性。故本文將選用SST (Menter)k-ω模型用于棒束通道的數值模擬計算。

1.2 一維User Code耦合方法

本文采用自行編寫的一維用戶程序(User Code)來模擬除棒束通道外的回路部分，包括預熱器、冷卻器與連接管路，如圖1所示。首先編寫一維用戶自定義程序，得到穩(wěn)態(tài)條件下系統(tǒng)回路的自然循環(huán)流量，和該自然循環(huán)流量對應的速度，進而將其設為CFD數值模擬的入口速度，最后模擬得到自然循環(huán)條件下棒束通道內的流動特性。

圖4 CFD不同模型與實驗數據壓降對比Fig.4 Comparison of pressure drop between different CFD models and experimental data

因此，編寫一維程序的目的就是找到實驗回路達到穩(wěn)定狀態(tài)的自然循環(huán)流量值，即驅動壓頭與總阻力壓頭相等時的流量值。設冷卻器與棒束通道中心高度差為L，棒束通道出口至冷卻器進口冷卻劑的密度為ρh，冷卻器出口到預熱器進口冷卻劑的密度為ρc，采用換熱中心假設，則驅動壓頭為ΔPd=(ρc-ρh)gL。總阻力壓頭ΔPz包括沿程阻力與局部阻力。沿程阻力即為各個管路的沿程阻力，具體計算方法為:

(1)

式中：λ為沿程阻力系數；l為管路長度；d為當量直徑；ρ為流過該管路的冷卻劑密度；A為管路流通面積；q為質量流量。其中沿程阻力系數λ的計算方法[12]為：

(2)

式中：a=0.094k0.225+0.53k;b=88k0.44;c=1.62k0.134;k=ε/d是相對粗糙度。

局部阻力包括棒束通道、預熱器、冷卻器處的局部阻力，以及各管路之間連接處，管路與棒束通道、預熱器、冷卻器連接處的局部阻力。局部阻力中除卻棒束通道內的阻力外，其他部分相對較小，可以忽略不計。那么，總阻力壓頭表示為ΔPz=ΔPf+ΔP-ρgH，ΔP為棒束通道進出口壓降，ρ為棒束通道的定性密度，H=800 mm為棒束通道長度。

因此，當實驗回路中驅動壓頭與總阻力壓頭相等時，即ΔPd=ΔPz時即可輸出自然循環(huán)回路質量流量q。具體程序框圖如圖5所示。

圖5 一維用戶程序框圖Fig.5 One-dimensional user code program chart

值得注意的是，將STAR-CCM+模擬的迭代次數輸入到一維程序中是為了編寫實現(xiàn)STAR-CCM+迭代次數超過一定步數(例如100步)后再進行上述循環(huán)過程的語句，目的是讓CFD達到相對收斂后再進行自然循環(huán)計算，從而避免計算結果的發(fā)散。

1.2.1 用戶庫

從上文可知，一維用戶程序需要將STAR-CCM+中的計算結果作為已知量，并且需要將用戶程序計算得到的自然循環(huán)速度傳遞給STAR-CCM+，使之作為棒束通道入口速度值。一維用戶程序與STAR-CCM+之間的數據傳遞需要通過事先把用戶自編程序編譯成鏈接庫，然后作為用戶庫加載到STAR-CCM+中的方式來實現(xiàn)。

1)用戶函數接口。

本文采用C++來編寫用戶函數，以實現(xiàn)得到自然循環(huán)回路流量的目的。通過庫注冊函數(uclib.cpp)指明用戶函數類型為標量場函數“ScalarFieldFunction”，以此在STAR-CCM+中棒束通道入口處通過設置場函數的方法來設置入口速度值。注冊方法可以通過調用ucfunc來實現(xiàn)，具體代碼為:

ucfunc(main, “ScalarFieldFunction”, “IterationVelocity”);

其中，main為用戶函數主函數名稱，IterationVelocity為顯示在STAR-CCM+下拉列表中的函數名。

另外，可通過調用ucarg注冊用戶函數所需的來自STAR-CCM+的參數，包括棒束通道進出口溫度、壓降，入口速度，迭代次數，具體代碼為:

ucarg(main, “Cell”, "$tinReport", sizeof(CoordReal));/*棒束通道入口溫度,K*/

ucarg(main, “Cell”, "$toutReport", sizeof(CoordReal));/*棒束通道出口溫度,K*/

ucarg(main, “Cell”, "$pjReport", sizeof(CoordReal));/*棒束通道壓降,Pa*/

ucarg(main, “Cell”, "$vReport", sizeof(CoordReal));/*棒束通道入口速度,m/s*/

ucarg(main, “Cell”, "$Iteration", sizeof(CoordReal));/*迭代次數*/

其中，Cell表示對于網格單元場的參數類型；Iteration為STAR-CCM+中原有的場函數，表示為迭代次數；tin、tout、pj、v為在STAR-CCM+中生成的報告，后綴加上Report可以同樣實現(xiàn)場函數的功能；size為變量組分表中元素所需的儲存(以字節(jié)為單位)，此大小可用于確保用戶函數的精度與STAR-CCM+的精度相匹配，所有的場函數參數都應設置為CoordReal，即表示為雙精度型浮點數據，另外盡管迭代次數Iteration為整型數據，也同樣需要設定為雙精度型浮點數據。值得注意的是，必須按照上述變量在用戶函數中的所需順序調用ucarg。本文中主函數的形參列表為：

void main(CoordReal *result, int size, CoordReal *tin, CoordReal *tout, CoordReal *pj, CoordReal *v, CoordReal *Iter)

其中，result為用戶函數的返回值的組分表，對于本文來說即為自然循環(huán)速度值;size為result組分表中的元素數量。

至此可以得到庫注冊函數uclib.cpp代碼如下：

#include "uclib.h"

void uclib()

{/*這里為上文所述的注冊用戶函數*/}

在與uclib.cpp相同的目錄中創(chuàng)建文件uclib.h，以聲明uclib.cpp中使用的函數。uclib.h文件中定義了在用戶函數中使用的變量和函數類型，它對于所有代碼都一樣，具體內容參見STAR-CCM+用戶指南[13]。

2)創(chuàng)建用戶庫。

至此，用戶庫源碼已生成完畢，包括uclib.h、 uclib.cpp、main.cpp, 由于在一維程序編寫過程中需要查物性參數，因此還需調用物性庫函數wasp97.h、wasp97.cpp。將上述文件同STAR-CCM+安裝目錄中的UserFunctions.lib編譯鏈接成為動態(tài)鏈接庫即可應用在STAR-CCM+中實現(xiàn)一維程序與三維CFD的耦合計算。

下面說明編譯方法。本文在Microsoft Visual C++ 2013上進行編譯，且Windows僅支持64位版本。打開VS2013 x64本機工具命令提示，將工作目錄定位到當前工作目錄，并使用以下命令將源程序uclib.cpp、main.cpp、wasp97.cpp編譯到對象文件中：

cl /MD /D_WINDOWS /DDOUBLE_PRECISION -c ***.cpp

即可將源代碼編譯成二進制目標文件uclib.obj, main.obj, wasp97.obj。

下面說明鏈接方法。在VS2013 x64本機工具命令提示中輸入如下命令：

link -dll /out:IterationVelocity.dll uclib.obj main.obj wasp97.obj "F:Program FilesCD-adapco13.04.011-R8STAR-CCM+13.04.011-R8starlibwin64intel16.3-r8libUserFunctions.lib"

鏈接成功后可在當前工作目錄中得到動態(tài)鏈接庫IterationVelocity.dll。將該用戶庫在STAR-CCM+模擬中加載，可發(fā)現(xiàn)場函數列表中新增場函數User IterationVelocity，將其設置為入口速度幅值即可實現(xiàn)一維User Code與三維CFD的耦合計算，研究自然循環(huán)條件下棒束通道的流動特性。

2 結果分析

通過上述方法可以對自然循環(huán)條件下的棒束通道進行數值模擬。當系統(tǒng)壓力為0.3 MPa時，在實驗操作中，可以通過提高預熱器加熱功率或棒束通道實驗段加熱功率來增加實驗回路的自然循環(huán)流量。對于數值模擬而言，為得到與實驗相對應的多組工況，提高預熱器加熱功率轉化為增加棒束通道實驗段入口溫度，提高棒束通道實驗段加熱功率可以通過直接在模擬中設置棒束熱流量來實現(xiàn)。

若保持棒束通道加熱功率不變，逐漸提高預熱器的加熱功率，可以得到系統(tǒng)回路自然循環(huán)流量隨棒束通道入口溫度變化關系。為方便與數值模擬結果對比，圖6展示了回路自然循環(huán)速度與棒束通道入口溫度關系曲線，其中模擬速度值為棒束通道進出口速度平均值。從圖中可以看出采用三維CFD軟件STAR-CCM+和一維自定義用戶程序耦合模擬可以很好地預測棒束通道內的自然循環(huán)速度值，可以認為STAR-CCM+與一維User Code耦合能夠研究棒束通道內的自然循環(huán)流動特性。針對圖1所示實驗過程中得到引壓管2、3間的壓降值ΔP2，同樣在模擬中得到相同部位的壓降，并與實驗值ΔP2進行對比，如圖6所示，可以看出二者符合較好，最大誤差在0.5%以內。

圖6 變預熱器功率實驗值與模擬值對比Fig.6 Comparison of experimental and simulation results of variable preheater power

同樣，若保持預熱器加熱功率不變，逐漸提高棒束通道的加熱功率，可以得到系統(tǒng)回路自然循環(huán)速度隨棒束通道熱流量的變化關系，如圖7所示。對于壓降ΔP2也采用與圖6同樣的處理方法，對比結果如圖7所示，發(fā)現(xiàn)二者仍然符合較好，最大誤差在0.5%以內?？梢钥闯鰺o論是改變預熱器功率或改變棒束通道加熱功率，STAR-CCM+與一維用戶程序耦合都能很好地預測棒束通道內的自然循環(huán)流動特性。

圖7 變棒束通道功率實驗值與模擬值對比Fig.7 Comparison of experimental and simulation results of variable rod bundle channel power

3 結論

1)在STAR-CCM+中對3×3棒束通道使用SST (Menter) K-Omega湍流模型進行數值模擬，自行編寫C++一維用戶程序求解系統(tǒng)回路其余部分的流動特性，再通過STAR-CCM+的用戶接口實現(xiàn)了二者的耦合。

2)在改變預熱器加熱功率和改變棒束通道加熱功率2種情況下將數值模擬計算得到的自然循環(huán)冷卻劑速度值和棒束通道內壓降結果與相應的實驗值進行對比，發(fā)現(xiàn)二者符合較好，可以說明STAR-CCM+與一維用戶程序耦合能夠很好的預測棒束通道的自然循環(huán)流動特性。

本文對一維用戶程序與三維STAR-CCM+之間的數據交互方法進行了說明，該耦合方法和實現(xiàn)流程為后續(xù)研究奠定了基礎。