杜為安,張 琦,楊繼行,閆超星,張 翼,劉成洋
(1.中國(guó)艦船研究設(shè)計(jì)中心,湖北 武漢 430061;2.上海交通大學(xué) 核科學(xué)與工程學(xué)院,上海 200240)
典型壓水反應(yīng)堆在一回路冷卻劑內(nèi)添加硼酸來補(bǔ)償剩余反應(yīng)性,硼酸溶液在反應(yīng)堆壓力容器以及一回路管道內(nèi)的分布特性直接影響堆芯功率[1-3]。當(dāng)堆芯供冷不足、堆芯內(nèi)產(chǎn)生大量水蒸氣時(shí),由于硼酸在蒸汽中溶解度較小,幾乎不含硼的水蒸氣經(jīng)過蒸汽發(fā)生器后將形成不含硼的冷凝液(稱為稀釋水團(tuán))。稀釋水團(tuán)將聚集在壓力容器入口上游的相關(guān)管路中[4],這部分稀釋水團(tuán)進(jìn)入壓力容器后將會(huì)和含硼冷卻劑混合,并最終進(jìn)入堆芯。堆芯應(yīng)急冷卻系統(tǒng)啟動(dòng)后,安全注射水箱中的高濃度硼溶液被注入壓力容器內(nèi)與一回路冷卻劑混合后進(jìn)入堆芯。安注冷卻劑和一回路冷卻劑之間還可能導(dǎo)致熱分層[5-6]。上述含硼溶液的混合過程直接關(guān)系到堆芯反應(yīng)性的變化,攪混不均勻可能導(dǎo)致堆芯反應(yīng)性激增或功率分布偏移,進(jìn)而威脅堆芯安全。
壓力容器的環(huán)形下降段位于堆芯入口上游,研究含硼溶液在環(huán)形下降段內(nèi)的混合過程對(duì)反應(yīng)堆安全具有重要意義[1,7]。國(guó)內(nèi)外已開展過含硼溶液混合相關(guān)的實(shí)驗(yàn)[2,8-9],受測(cè)量技術(shù)限制,現(xiàn)階段的相關(guān)實(shí)驗(yàn)只能獲得環(huán)形下降段內(nèi)局部位置的平均硼濃度分布。計(jì)算流體力學(xué)(CFD)軟件成本較低、可重復(fù)性高,并可對(duì)三維場(chǎng)的流動(dòng)與混合過程進(jìn)行分析,被廣泛應(yīng)用于反應(yīng)堆熱工水力的相關(guān)研究。本文基于計(jì)算流體力學(xué)軟件ANSYS Fluent 18.0對(duì)環(huán)形下降段內(nèi)的流動(dòng)與混合特性進(jìn)行分析。
本文基于不可壓縮流體假設(shè),忽略溫度變化對(duì)濃度擴(kuò)散的影響,控制方程如下。
1) 連續(xù)方程
連續(xù)方程即質(zhì)量守恒方程,其微分形式為:
(1)
式中:u、v、w為x、y、z方向速度分量,m/s;ρ為密度,kg/m3;t為時(shí)間,s。
2) 動(dòng)量方程
動(dòng)量方程實(shí)際上是牛頓第二定律,x、y、z方向的動(dòng)量方程為:
(2)
(3)
(4)
式中:p為流體微元體上的壓力,Pa;τxx、τxy、τxz等為黏性應(yīng)力τ的分量,Pa;fx、fy、fz為x、y、z方向的單位質(zhì)量力,m/s2,若質(zhì)量力只有重力,z軸豎直向上,則fx=fy=0,fz=-g。
3) 能量方程
包含有熱交換的系統(tǒng)必須滿足能量守恒定律,其本質(zhì)為熱力學(xué)第一定律,表達(dá)式為:
(5)
式中:E為流體微團(tuán)總焓值,J/kg,由內(nèi)能、勢(shì)能和動(dòng)能組成,E=h-p/ρ+u2/2;hj為組分j的焓值,J/kg;Jj為組分j的擴(kuò)散通量;k為流體熱傳導(dǎo)系數(shù),W/(m·K);Sh為包括化學(xué)反應(yīng)及其他體積熱源項(xiàng)。
4) 組分質(zhì)量守恒方程
在存在質(zhì)的交換或多種化學(xué)組分的特性控制系統(tǒng)中,每一組分都要遵循組分質(zhì)量守恒定律。對(duì)于一特定的系統(tǒng),組分質(zhì)量守恒定律可描述為:系統(tǒng)內(nèi)某組分對(duì)時(shí)間的變化率等于通過系統(tǒng)界面的凈擴(kuò)散流量和化學(xué)反應(yīng)生成率之和。
由組分守恒定律,得到組分s的組分質(zhì)量守恒方程:
(6)
式中:cs、ρcs為組分s的體積濃度和質(zhì)量濃度;Ds為組分s的擴(kuò)散系數(shù);Ss為源項(xiàng)。式(6)從左到右分別為時(shí)間變化率、對(duì)流項(xiàng)、擴(kuò)散項(xiàng)和源項(xiàng)。
本文參照哈爾濱工程大學(xué)激光診斷實(shí)驗(yàn)室的硼濃度混合分析實(shí)驗(yàn)[10]建立相應(yīng)的物理模型,以便于對(duì)比實(shí)驗(yàn)結(jié)果與數(shù)值模擬結(jié)果。表1為壓力容器物理模型的相關(guān)參數(shù)。圖1為壓力容器物性模型的俯視圖和正視圖。流體經(jīng)入口管進(jìn)入環(huán)形下降段,流入下腔室后從下腔室上部的7×7矩形通道流出。下腔室并不是本文關(guān)注的重點(diǎn),在保證流體流動(dòng)方向不變的情況下,簡(jiǎn)化下腔室內(nèi)構(gòu)建,僅保留半球形結(jié)構(gòu)和上部孔板。
表1 壓力容器模型參數(shù)Table 1 Parameter of reactor model
圖1 壓力容器簡(jiǎn)化模型Fig.1 Mock RPV model
采用ANSYS ICEM對(duì)模型進(jìn)行網(wǎng)格劃分,網(wǎng)格結(jié)構(gòu)如圖2所示。將壓力容器模型分為入口管段、環(huán)形下降段、下腔室和矩形出口通道4個(gè)主要區(qū)域,各區(qū)域均采用六面體結(jié)構(gòu)化網(wǎng)格,各區(qū)域交接面的節(jié)點(diǎn)一一對(duì)應(yīng)。經(jīng)無關(guān)性驗(yàn)證后(圖3),最終網(wǎng)格數(shù)量約為127萬,網(wǎng)格質(zhì)量為0.55。
Rohde等[11]結(jié)合硼混合實(shí)驗(yàn)數(shù)據(jù)進(jìn)行湍流模型驗(yàn)證,并指出不考慮浮升力影響的情況下,SST模型可準(zhǔn)確計(jì)算壓力容器內(nèi)硼濃度分布。SST模型對(duì)壁面附近的低雷諾數(shù)流動(dòng)區(qū)域采用標(biāo)準(zhǔn)k-ω模型,對(duì)充分發(fā)展的自由流動(dòng)區(qū)域采用標(biāo)準(zhǔn)k-ε模型,因此很大程度上降低了標(biāo)準(zhǔn)k-ω模型對(duì)自由剪切流的敏感程度[12-13]。SST模型的輸運(yùn)方程如下:
圖2 壓力容器模型網(wǎng)格劃分Fig.2 Meshing of RPV model
圖3 網(wǎng)格無關(guān)性驗(yàn)證Fig.3 Independence validation of mesh number
(7)
式中:Gk為由時(shí)均速度梯度引起的湍流動(dòng)能產(chǎn)生項(xiàng);Gω為湍流渦頻率產(chǎn)生項(xiàng);Yk和Yω分別為湍流動(dòng)能及湍流渦頻率的耗散項(xiàng);Dω為交叉擴(kuò)散項(xiàng)。
為驗(yàn)證CFD計(jì)算結(jié)果的準(zhǔn)確性,將數(shù)值模擬和實(shí)驗(yàn)得到的入口橫截面平均質(zhì)量分?jǐn)?shù)進(jìn)行對(duì)比,結(jié)果如圖4所示??煽闯?,入口橫截面的平均質(zhì)量分?jǐn)?shù)隨注射時(shí)間的增加而增大,在10 s后增長(zhǎng)速度明顯減緩。CFD計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果較為一致,平均相對(duì)誤差為3.41%。說明本文計(jì)算結(jié)果具有較高的準(zhǔn)確性。
圖4 數(shù)值模擬結(jié)果與實(shí)驗(yàn)對(duì)比Fig.4 Comparison between CFD simulation and experiment
開啟混合物模塊,自定義硼酸溶液和水溶液的混合物,通過編寫UDF設(shè)置二者的組分占比,入口流速分別設(shè)置為0.2、0.3、0.4 m/s,硼酸質(zhì)量分?jǐn)?shù)c設(shè)置為:t≤0 s時(shí),c=0;t>0 s時(shí),c=0.1。系統(tǒng)運(yùn)行壓力為0.1 MPa,混合物的物性設(shè)置為常壓下水的物性。忽略硼示蹤劑對(duì)流體的物性影響,流體密度為998.2 kg/m3,運(yùn)動(dòng)黏度為0.001 Pa·s。入口邊界條件設(shè)置為速度入口,出口邊界條件設(shè)置為壓力出口,壁面類型選擇無滑移壁面邊界條件。模擬混合過程前首先將壓力容器內(nèi)的全部流體定義為水,t>0 s時(shí)含水混合物自入口注入壓力容器。
1) 二維速度場(chǎng)分析
圖5 水平方向速度矢量分布Fig.5 Transverse velocity vector
圖5為橫截面的速度矢量分布,對(duì)應(yīng)的入口雷諾數(shù)為11 448??煽闯?,在0 cm橫截面周向位置[10°,30°]和[-30°,-10°]區(qū)域存在較大尺度的渦結(jié)構(gòu),由于入口流體與內(nèi)壁發(fā)生劇烈碰撞,流體呈放射狀散開,在靠近入口管道的區(qū)域形成回流并與環(huán)腔內(nèi)流體混合。Eltayeb等[14]通過可視化實(shí)驗(yàn)在該區(qū)域觀測(cè)到了旋渦的形成與瓦解過程,旋渦的尺寸與位置存在一定的隨機(jī)性。而數(shù)值模擬結(jié)果僅能觀測(cè)到1個(gè)尺度較大的渦,這是由于兩方程渦黏模型基于平均化假設(shè),在瞬態(tài)條件下無法計(jì)算出脈動(dòng)值[15]。在單回路流動(dòng)實(shí)驗(yàn)的初始階段(0~2 s),示蹤劑在上述區(qū)域的擴(kuò)散過程與速度矢量分布相對(duì)應(yīng),染色劑與內(nèi)壁發(fā)生碰撞后形成逆流向外壁發(fā)展。隨后流體沿環(huán)形區(qū)域周向位置流動(dòng),期間不斷與環(huán)腔壁面發(fā)生碰撞,水平方向的速度隨周向位置的延伸而衰減,在[90°,130°]和[-130°,-90°]的區(qū)域最為明顯。在水平位置較低的-15 cm橫截面,[45°,135°]和[-45°,-135°]區(qū)域的流體流動(dòng)不受上游渦結(jié)構(gòu)的影響,速度分布較為平穩(wěn)。
圖6、7分別為環(huán)形下降段內(nèi)水平方向和豎直方向的速度分布。如前文所述,流體進(jìn)入環(huán)腔后沿周向位置擴(kuò)散,期間不斷與壁面發(fā)生碰撞,水平方向動(dòng)能不斷損失,水平方向的速度隨周向位置的增加而衰減。注意到-10~-25 cm截面水平方向速度先增加后減小,其中以圖6a最為明顯,這是由于靠近入口管道的區(qū)域受放射狀回流的影響。部分入口流體與內(nèi)壁碰撞后將沿豎直方向進(jìn)入環(huán)腔下部空間;環(huán)腔上部空間的入口流體受質(zhì)量力影響也會(huì)逐漸沿豎直方向沉降,可看出水平位置越低的截面,豎直方向的速度分量越大。
2) 三維流線圖分析
圖8為單回路穩(wěn)態(tài)流動(dòng)條件下環(huán)形下降段和下腔室內(nèi)的流線圖。流體經(jīng)入口管進(jìn)入壓力容器后與環(huán)腔內(nèi)壁發(fā)生劇烈碰撞并形成回流,部分回流受邊界條件的束縛沿內(nèi)壁呈放射狀散開。隨后入口流體分成兩股沿著環(huán)形下降段周向位置流動(dòng),豎直方向的速度分量隨周向位置的延伸而增加,到達(dá)入口對(duì)側(cè)后環(huán)腔左右兩股流體開始匯合。
圖6 水平方向速度分布Fig.6 Transverse velocity distribution
圖7 豎直方向速度分布Fig.7 Longitudinal velocity distribution
圖8 環(huán)腔與下腔室內(nèi)流線圖Fig.8 Streamlines in down-comer and lower plenum
流線圖可反映微元流動(dòng)介質(zhì)在流場(chǎng)中的軌跡。從圖8可看出,在入口正下方和環(huán)腔上部空間存在流動(dòng)停滯區(qū),該區(qū)域隨入口雷諾數(shù)的增加而收縮。在主泵恢復(fù)工作的工況下,主泵轉(zhuǎn)速由控制器驅(qū)動(dòng)呈點(diǎn)動(dòng)式增加而非線性增加,持續(xù)的低速轉(zhuǎn)動(dòng)可能導(dǎo)致稀釋水團(tuán)滯留在環(huán)腔上部空間,與冷卻劑交混過程緩慢。
1) 橫向濃度場(chǎng)分布
將數(shù)值模擬得到的橫截面混合物質(zhì)量分?jǐn)?shù)分布與文獻(xiàn)[10]所報(bào)道的示蹤劑濃度分布進(jìn)行比較,結(jié)果如圖9、10所示。從數(shù)值模擬結(jié)果可看出,混合流體進(jìn)入環(huán)腔后與內(nèi)壁發(fā)生碰撞,隨后沿著周向位置發(fā)展并不斷與環(huán)腔內(nèi)流體混合;周向位置越大,混合物擴(kuò)散的速度越慢。在圖9、10中,數(shù)值模擬圖像中的比色卡為混合物份額,而實(shí)驗(yàn)結(jié)果圖像的比色卡為染色劑(RhB)的濃度,二者的擴(kuò)散趨勢(shì)高度一致,進(jìn)一步驗(yàn)證了模型的可靠性。SST模型能得到濃度擴(kuò)散整體趨勢(shì),但無法獲得湍流擴(kuò)散的細(xì)節(jié)。在數(shù)值模擬結(jié)果中可觀測(cè)到入口位置內(nèi)壁和外壁之間存在一定的濃度梯度,而實(shí)驗(yàn)結(jié)果顯示該區(qū)域的混合較為充分,濃度梯度較小。
a——數(shù)值模擬質(zhì)量分?jǐn)?shù)分布;b——LIF實(shí)驗(yàn)染色劑濃度分布圖9 0 cm橫截面實(shí)驗(yàn)與數(shù)值模擬結(jié)果對(duì)比Fig.9 Contrast between numerical simulation and experimental results at 0 cm cross section
a——數(shù)值模擬質(zhì)量分?jǐn)?shù)分布;b——LIF實(shí)驗(yàn)染色劑濃度分布圖10 -10 cm橫截面實(shí)驗(yàn)與數(shù)值模擬結(jié)果對(duì)比Fig.10 Contrast between numerical simulation and experimental results at -10 cm cross section
圖11為周向位置質(zhì)量分?jǐn)?shù)數(shù)值模擬與濃度實(shí)驗(yàn)結(jié)果。從圖11a可觀測(cè)到,混合物進(jìn)入環(huán)腔后在入口位置持續(xù)出現(xiàn)1個(gè)質(zhì)量分?jǐn)?shù)的峰值,由于入口流體在水平方向通過的距離較短,且尚未與環(huán)形下降段內(nèi)外壁發(fā)生碰撞,混合特性較差。隨后高濃度硼溶液沿周向位置流動(dòng)并不斷與水發(fā)生混合,環(huán)腔內(nèi)混合物所占區(qū)域不斷增加,如圖11b所示??勺⒁獾剑趖>5 s后,周向角為[30°,110°]和[-30°,-110°]區(qū)域的質(zhì)量分?jǐn)?shù)平穩(wěn)增長(zhǎng),說明該區(qū)域的混合物已與水充分混合。從圖11可看出,周向位置質(zhì)量分?jǐn)?shù)不斷展平,質(zhì)量分布區(qū)域均勻。由于壓力容器模型加工精度所限環(huán)腔兩側(cè)濃度分布輕微不對(duì)稱,整體上,數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果較為一致,實(shí)際流動(dòng)過程染色劑發(fā)展的邊緣位置存在一定脈動(dòng)。
2) 三維濃度場(chǎng)分布
圖12為硼酸溶液注入環(huán)形下降段后的質(zhì)量分?jǐn)?shù)分布。含硼溶液進(jìn)入壓力容器后呈放射狀散開,并繞環(huán)腔壁邊界流動(dòng);混合物到達(dá)周向位置較大的區(qū)域后水平方向的速度衰減,其輸運(yùn)特性過渡為沿豎直方向的沉降,以圖12b中t=24 s時(shí)刻最為明顯。從圖12可看出,單回路流動(dòng)混合時(shí),入口正下方存在1個(gè)流動(dòng)停滯區(qū),在整個(gè)混合過程中模擬硼酸溶液無法進(jìn)入到流動(dòng)停滯區(qū)??梢姡葧r(shí)間下,入口混合流體雷諾數(shù)越大,混合物所分布的區(qū)域相應(yīng)越大,整個(gè)環(huán)腔內(nèi)的流通停滯區(qū)相應(yīng)收縮。文獻(xiàn)[10]所報(bào)道的示蹤劑濃度擴(kuò)散特性與該現(xiàn)象一致。
a——周向位置質(zhì)量分?jǐn)?shù)分布數(shù)值模擬;b——周向位置濃度分布LIF實(shí)驗(yàn)結(jié)果圖11 數(shù)值模擬與實(shí)驗(yàn)結(jié)果定量對(duì)比Fig.11 Quantitative contrast between numerical simulation and experimental results
a——Re=5 833;b——Re=11 448圖12 壓力容器三維質(zhì)量分?jǐn)?shù)分布Fig.12 Three-dimensional mass fraction distribution in RPV
本文基于SST模型對(duì)環(huán)形下降段內(nèi)的硼混合問題進(jìn)行了分析,得到如下結(jié)論。
1) 單回路條件下,環(huán)形下降段上部空間水平方向的速度梯度較大,下部空間流體流動(dòng)較為穩(wěn)定;環(huán)形下降段內(nèi)存在流動(dòng)停滯區(qū),且該區(qū)域隨入口雷諾數(shù)的增加而收縮。
2) 含硼溶液從單一入口進(jìn)入壓力容器后將沿環(huán)形下降段周向輸運(yùn),期間水平方向動(dòng)能不斷損失,在環(huán)腔內(nèi)的輸運(yùn)特性逐漸過渡為沿豎直方向的沉降。