STAR-CCM+與一維用戶程序耦合方法

張銀星， 高璞珍， 何曉強(qiáng)， 陳沖， 林宇琦， 劉怡雯

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

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

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

1 三維與一維耦合方法

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

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

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

1.1 三維CFD數(shù)值模擬

棒束通道實(shí)驗(yàn)段結(jié)構(gòu)如圖2所示。棒束通道實(shí)驗(yàn)段長度為800 mm，利用STAR-CCM+內(nèi)的3D-CAD模型可直接進(jìn)行建模。

1.1.1 網(wǎng)格獨(dú)立性分析

本文采用多面體網(wǎng)格實(shí)現(xiàn)棒束通道主體結(jié)構(gòu)的網(wǎng)格劃分，采用棱柱層網(wǎng)格模型應(yīng)用在所有的壁面邊界。圖3給出棒束通道內(nèi)生成的三維網(wǎng)格。

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

圖2 棒束通道結(jié)構(gòu)Fig.2 Rod bundle channel structure

圖3 棒束通道網(wǎng)格示意Fig.3 Rod bundle channel grid diagram

表1 網(wǎng)格無關(guān)性驗(yàn)證Table 1 Grid independence verification

1.1.2 模型驗(yàn)證

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

1.2 一維User Code耦合方法

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

圖4 CFD不同模型與實(shí)驗(yàn)數(shù)據(jù)壓降對(duì)比Fig.4 Comparison of pressure drop between different CFD models and experimental data

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

(1)

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

(2)

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

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

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

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

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

1.2.1 用戶庫

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

1)用戶函數(shù)接口。

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

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

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

另外，可通過調(diào)用ucarg注冊(cè)用戶函數(shù)所需的來自STAR-CCM+的參數(shù)，包括棒束通道進(jìn)出口溫度、壓降，入口速度，迭代次數(shù)，具體代碼為:

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));/*迭代次數(shù)*/

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

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

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

至此可以得到庫注冊(cè)函數(shù)uclib.cpp代碼如下：

#include "uclib.h"

void uclib()

{/*這里為上文所述的注冊(cè)用戶函數(shù)*/}

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

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

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

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

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

即可將源代碼編譯成二進(jìn)制目標(biāo)文件uclib.obj, main.obj, wasp97.obj。

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

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"

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

2 結(jié)果分析

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

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

圖6 變預(yù)熱器功率實(shí)驗(yàn)值與模擬值對(duì)比Fig.6 Comparison of experimental and simulation results of variable preheater power

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

圖7 變棒束通道功率實(shí)驗(yàn)值與模擬值對(duì)比Fig.7 Comparison of experimental and simulation results of variable rod bundle channel power

3 結(jié)論

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

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

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