迷宮型節(jié)流件壓降的數(shù)值模擬研究

宋 怡，陸道綱，秦亥琦，鐘達(dá)文，曹 瓊

(華北電力大學(xué) 非能動核能安全技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 102206)

池式鈉冷快堆主冷卻系統(tǒng)的流道包括主循環(huán)流道和五個(gè)分支流道。在分支流道中，迷宮型節(jié)流件被安裝在一些重要的部件前后來調(diào)節(jié)壓降和控制流量，例如在主循環(huán)泵的旁通流量計(jì)中，鈉從泵后側(cè)引出的接管流出，依次流經(jīng)迷宮型節(jié)流件、連接管道、流量計(jì)、出口管，最終流進(jìn)主循環(huán)泵的泵箱。除此之外，泵支承冷卻系統(tǒng)、堆容器冷卻系統(tǒng)、堆內(nèi)電離室冷卻系統(tǒng)中也用到了迷宮型節(jié)流件。當(dāng)冷卻劑鈉流過這些迷宮型節(jié)流件時(shí)，鈉溫均為360 ℃[1]。由于在設(shè)計(jì)和運(yùn)行過程中反應(yīng)堆的熱工水力參數(shù)經(jīng)常被調(diào)整，因此通過這些迷宮型節(jié)流件的流量和需要實(shí)現(xiàn)的壓降也會有所調(diào)整。為了更加合理地挑選迷宮型節(jié)流件，需要較為準(zhǔn)確地評估它的壓降。

所有三維迷宮型節(jié)流件均由厚度相同、高度相同且高于管道半徑的節(jié)流片組成，相鄰節(jié)流片間距相同且呈180°交錯(cuò)排列。節(jié)流件放置在圓形管道中，示意圖如圖1所示。

圖1 迷宮型節(jié)流件示意圖

評估壓降通常采用實(shí)驗(yàn)的方法，但是實(shí)驗(yàn)不僅費(fèi)時(shí)費(fèi)力，而且難以觀察流動細(xì)節(jié)，近年來，數(shù)值模擬已經(jīng)成為計(jì)算壓降的常用方法。國內(nèi)外許多學(xué)者用數(shù)值模擬方法研究了各種節(jié)流件的流動和阻力特性。Abdulrazaq A.使用realizablek-ε模型研究了兩個(gè)相似的連續(xù)排列的斜邊孔板的流動特性。發(fā)現(xiàn)收縮位置與第一個(gè)孔有關(guān)，與第二個(gè)孔無關(guān)，且隨孔徑比的變化而變化[2]。Ehsan S應(yīng)用數(shù)值模擬方法，優(yōu)化渦流管的參數(shù)和幾何形狀，獲得了能夠?qū)崿F(xiàn)最高制冷效率的節(jié)流件錐長[3]。Guohui Gan采用CFD方法預(yù)測了方形薄板上單孔和多孔的壓力損失特性，結(jié)論是多孔的壓力損失系數(shù)比單孔高約0.5[4]。Ivan Okhotnikov對一種新型單級節(jié)流閥進(jìn)行了CFD數(shù)值模擬，預(yù)測了由于穩(wěn)定流動引起的扭矩和壓降，分析了開口面積、流量系數(shù)和射流角度的圖像，證明了該設(shè)計(jì)在高流量的流體動力系統(tǒng)中使用的有效性[5]。王海民通過數(shù)值模擬方法研究了由一系列直角轉(zhuǎn)彎組成的曲折迷宮通道中的流速和壓力分布，證明了由串聯(lián)通道產(chǎn)生的多級壓降可以在控制流速合理的同時(shí)防止空化現(xiàn)象的產(chǎn)生[6]。Han Xu模擬了帶有環(huán)形塞的止回閥內(nèi)的湍流流動，通過與實(shí)驗(yàn)結(jié)果對比，證明RNGk-ε模型足以預(yù)測該閥的水動力特性[7]。

盡管數(shù)值模擬很便捷，但它的準(zhǔn)確性受很多因素影響，尤其是湍流模型的選擇。目前最適合迷宮型節(jié)流件數(shù)值模擬的湍流模型尚不明確，數(shù)值模擬的結(jié)果也缺少實(shí)驗(yàn)驗(yàn)證。

另一方面，迷宮型節(jié)流件的壓降是由多種因素決定的，包括相鄰節(jié)流板的間距、節(jié)流板的高度、節(jié)流板的數(shù)目和流過節(jié)流件的流量。目前關(guān)于迷宮型節(jié)流件的壓降和這些因素之間的關(guān)系尚不明確。多元回歸分析方法是變量關(guān)系研究常用的方法，通過它擬合出顯式公式不僅可以較為直觀地表明變量之間的關(guān)系，而且可以通過已有數(shù)據(jù)來預(yù)測更多數(shù)據(jù)。Qi Sun使用計(jì)算流體力學(xué)方法，分析了顆粒雷諾數(shù)、流動角和阻力系數(shù)百分比對高溫氣冷堆中簡化的石墨粉塵顆粒阻力系數(shù)的影響。并基于Rosendahl提出的框架采用多元回歸方法擬合了更精確的阻力系數(shù)公式。他還應(yīng)用最小二乘法求出相關(guān)參數(shù)，顯著降低了預(yù)測阻力系數(shù)的誤差[8]。Karim Nassar對海堤的水動力性能進(jìn)行了實(shí)驗(yàn)，分析了各種因素對海堤水動力性能的影響，利用多元非線性回歸方法開發(fā)了新的經(jīng)驗(yàn)公式。不僅用于理解擬建海堤的水動力特性，并且預(yù)測了某些特性。基于實(shí)驗(yàn)結(jié)果，這些公式用軟件SPSS進(jìn)行了驗(yàn)證[9]。

基于上述認(rèn)識，本文將研究計(jì)算迷宮型節(jié)流件壓降的數(shù)值模擬方法，并且通過數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對比，選擇最可信的湍流模型。本文將使用經(jīng)過驗(yàn)證的數(shù)值模擬方法，分析迷宮型節(jié)流件的壓降隨各種因素變化的趨勢，最終通過多元回歸方法擬合壓降與這些因素之間關(guān)系的經(jīng)驗(yàn)公式。

1 迷宮型節(jié)流件的數(shù)值模擬方法及其驗(yàn)證

1.1 迷宮型節(jié)流件的數(shù)值模擬方法

將通過實(shí)驗(yàn)驗(yàn)證的三維迷宮型節(jié)流件，放置在內(nèi)徑123 mm的圓形管道中。它由兩塊節(jié)流片組成，每塊節(jié)流片厚度為7 mm，高度為78 mm，相鄰節(jié)流片間距為16 mm，呈180°交錯(cuò)排列，如圖2所示。為了使數(shù)值模擬的幾何模型與實(shí)驗(yàn)的幾何模型完全一致，我們將節(jié)流件的入口延長615 mm，出口延長1 230 mm。

迷宮型節(jié)流件的幾何模型在Solid Works中創(chuàng)建，然后將其導(dǎo)入ICEM CFD生成非結(jié)構(gòu)網(wǎng)格，并在近壁區(qū)生成良好的邊界層，生成的網(wǎng)格將導(dǎo)入Fluent進(jìn)行計(jì)算。Fluent的計(jì)算過程基于有限元方法，也就是將連續(xù)的流體域劃分為離散單元，再對每個(gè)單元求解方程。如果單元數(shù)太少，難以模擬流動細(xì)節(jié)，計(jì)算結(jié)果將產(chǎn)生較大誤差；如果單元數(shù)太多，計(jì)算機(jī)負(fù)擔(dān)過重，會影響計(jì)算速度。因此，在選擇湍流模型以前，我們分析了網(wǎng)格獨(dú)立性，確定了能夠保證計(jì)算精度的最小網(wǎng)格數(shù)。網(wǎng)格獨(dú)立性分析(使用標(biāo)準(zhǔn)k-ε模型，雷諾數(shù)為852 744)的結(jié)果如圖3所示，220萬單元數(shù)的網(wǎng)格被選為最佳網(wǎng)格。

圖2 數(shù)值模擬的幾何模型

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

由于目前最適合迷宮型節(jié)流件數(shù)值模擬的湍流模型尚不明確，因此我們運(yùn)用了多個(gè)湍流模型來重復(fù)數(shù)值模擬過程，包括標(biāo)準(zhǔn)k-ε模型、標(biāo)準(zhǔn)k-ω模型、可實(shí)現(xiàn)的k-ε模型,RNGk-ε模型,SSTk-ω模型和Spalart-Allmaras模型。它們的控制方程可以從已經(jīng)發(fā)表的文獻(xiàn)中查閱[10-18]。將這些湍流模型計(jì)算得到的數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果對比，得到最適合迷宮型節(jié)流件的湍流模型。

數(shù)值模擬中使用的流體是360 ℃的液態(tài)鈉，與快堆主冷卻系統(tǒng)支路中的鈉溫一致。液鈉被認(rèn)為是單相的、不可壓縮的，它的物理性質(zhì)不隨溫度變化。Fluent中的邊界條件如表1所示。

表1 邊界條件

1.2 驗(yàn)證

1.2.1 實(shí)驗(yàn)介紹

實(shí)驗(yàn)設(shè)備包括水箱、泵、閥門、加熱器、流量計(jì)、差壓計(jì)和節(jié)流件的實(shí)驗(yàn)段，如圖4所示。其中水箱用來儲存純凈水，并為實(shí)驗(yàn)環(huán)路提供水源。加熱器用來維持環(huán)路中的水溫。

圖4 實(shí)驗(yàn)裝置和環(huán)路示意圖

節(jié)流件的實(shí)驗(yàn)段用不銹鋼鍛造，并通過法蘭連接到實(shí)驗(yàn)環(huán)路。節(jié)流件的前后均布置有壓力引管，差壓計(jì)通過壓力引管連接到實(shí)驗(yàn)環(huán)路。為了確保差壓變送器所測的壓降不受流量分布不均的影響，我們將取壓點(diǎn)布置在相同高度，且采用環(huán)形四點(diǎn)取壓方式。實(shí)驗(yàn)中，流體從節(jié)流件底部的入口流入，流經(jīng)節(jié)流件后從頂部的出口流出，通過調(diào)節(jié)閥門的開度來改變流量。流量計(jì)和差壓計(jì)自動收集體積流量和壓降的數(shù)據(jù)，經(jīng)信號轉(zhuǎn)換模塊轉(zhuǎn)換后顯示并儲存在電腦端。

鑒于液態(tài)鈉價(jià)格高昂及其危險(xiǎn)性，我們使用水來代替鈉進(jìn)行實(shí)驗(yàn)。為了能夠?qū)?shí)驗(yàn)結(jié)果應(yīng)用于原始模型，實(shí)驗(yàn)需要滿足相似理論，包括幾何相似、運(yùn)動相似和動力相似。

采用84 ℃的水進(jìn)行實(shí)驗(yàn)，得到壓降隨體積流量變化的實(shí)驗(yàn)結(jié)果如圖5所示。

圖5 壓降隨體積流量變化的實(shí)驗(yàn)結(jié)果

1.2.2 數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對比

我們使用不同湍流模型重復(fù)進(jìn)行了迷宮型節(jié)流件的數(shù)值模擬過程，質(zhì)量流量范圍是11～31 kg/s，雷諾數(shù)范圍是453 235～1 199 368。數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對比如圖6所示。在整個(gè)流量范圍內(nèi)，使用標(biāo)準(zhǔn)k-ε模型得到的數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果的誤差均在3%以內(nèi)。因此，我們選擇標(biāo)準(zhǔn)k-ε模型為計(jì)算迷宮型節(jié)流件壓降最可信的湍流模型。

2 迷宮型節(jié)流件水力特性的分析

使用上述經(jīng)過驗(yàn)證的數(shù)值模擬方法，我們分析了迷宮型節(jié)流件的水力特性，節(jié)流件由兩到四塊節(jié)流片組成，放置在內(nèi)徑為100 mm的管道中。

圖6 數(shù)值模擬與實(shí)驗(yàn)得到的壓降對比

2.1 空化現(xiàn)象及發(fā)生條件

在鈉冷快堆中，一旦在任何部件內(nèi)出現(xiàn)氣泡，氣泡就會混雜在液鈉中，隨液鈉流向其他部件，這會惡化冷卻劑的傳熱性能，對反應(yīng)堆造成潛在的危害。當(dāng)液鈉的壓力低于該溫度下的飽和蒸氣壓時(shí)，就會出現(xiàn)空化現(xiàn)象并產(chǎn)生氣泡。在節(jié)流件中，由于局部壓力的快速下降，空化現(xiàn)象更容易發(fā)生。由于迷宮型節(jié)流件被安裝在鈉冷快堆的主冷卻系統(tǒng)中，因此流過它的液鈉的絕對壓力應(yīng)當(dāng)高于360 ℃下的飽和蒸氣壓，防止空化現(xiàn)象的發(fā)生。

液態(tài)鈉的飽和蒸氣壓與溫度的關(guān)系可以通過多種方法獲得。Golden和Tokar對這些方法進(jìn)行了分析討論[19]并獲得了溫度低于1 250 K時(shí)的最佳方程，該方程由Ditchburn和Gilmour給出[20]：

(1)

式中：Ps——飽和蒸氣壓,Pa;

T——溫度,K。

當(dāng)溫度為360 ℃時(shí)，根據(jù)式(1)可以計(jì)算出飽和蒸氣壓為1.36×10-4Pa。

在迷宮型節(jié)流件中，流體每流過一個(gè)節(jié)流片，壓力都會迅速下降然后緩慢回升。另外，除了最后一個(gè)節(jié)流片以外，每個(gè)節(jié)流片的后側(cè)都有一個(gè)小漩渦，漩渦中心的壓力低于周圍區(qū)域的壓力。因此壓力最小的位置位于倒數(shù)第二塊節(jié)流板的后側(cè)。(見圖7的A、B、C點(diǎn))

圖7 流線圖與絕對壓力云圖

2.2 壓降隨各因素變化的趨勢

迷宮型節(jié)流件的壓降是由很多因素決定的，包括相鄰節(jié)流板的間距(D)、節(jié)流板的高度(H)、節(jié)流板的數(shù)目(N)和流過節(jié)流件的流量(M)，如圖8所示。本文分別分析壓降隨這些因素變化的趨勢。

圖8 影響迷宮型節(jié)流件壓降的各因素

在迷宮型節(jié)流件中，H和D決定了流體的流通面積，因此我們首先分析壓降與H、D之間的關(guān)系。隨著H減小或D增大，流通面積都會增大，進(jìn)而壓降減小。當(dāng)流體面積大到一定程度再繼續(xù)增大時(shí)，壓降幾乎不變。另外，當(dāng)流通面積太小時(shí)，空化現(xiàn)象更易發(fā)生。

圖9展示了壓降隨D、H的變化趨勢。當(dāng)D小于引起空化的臨界間距時(shí)，液鈉的最小絕對壓力會低于該溫度下的飽和蒸氣壓，將導(dǎo)致空化并造成危險(xiǎn)；當(dāng)D小于25 mm，隨著D增加，壓降顯著減小；當(dāng)D大于25 mm時(shí)，隨著D增加，壓降幾乎不變。因此，D的最佳范圍是從引起空化的臨界間距(Dmin)到25 mm。當(dāng)H小于65 mm時(shí)，隨著H增大，壓降緩慢增加；當(dāng)H大于65 mm時(shí)，隨著H增大，壓降顯著增加。

圖9 壓降隨D、H的變化趨勢

圖10展示了壓降隨N、M的變化趨勢。由于迷宮型節(jié)流件的降壓形式是多級降壓，因此隨著N增加，壓降逐漸增大。隨著M增加，單位時(shí)間內(nèi)流過節(jié)流件的流體體積增大，因此壓降增大。

圖10 壓降隨N、M的變化趨勢

2.3 迷宮型節(jié)流件壓降的經(jīng)驗(yàn)公式

由于在設(shè)計(jì)和運(yùn)行過程中反應(yīng)堆的熱工水力參數(shù)經(jīng)常被調(diào)整，因此通過這些迷宮型節(jié)流件的流量和需要實(shí)現(xiàn)的壓降也會有所調(diào)整。我們定量分析了大量數(shù)值模擬結(jié)果，用多元非線性回歸分析方法擬合了經(jīng)驗(yàn)公式，擬合出的經(jīng)驗(yàn)公式可以被用作評估迷宮型節(jié)流件壓降的參考。經(jīng)驗(yàn)公式的適用范圍是管道內(nèi)徑100 mm，H為60～80 mm，D為Dmin—25 mm，N為2～4，M為4～12 kg/s。

由于四個(gè)自變量(N,M,H,D)在物理意義上相互獨(dú)立，因此我們分別畫出了表示壓降P與這些自變量之間關(guān)系的散點(diǎn)圖，如圖11所示。

圖11 壓降與自變量關(guān)系的散點(diǎn)圖

如圖11所示，壓降與N、M呈現(xiàn)一定的線性關(guān)系；此外，由于N和M維度的樣本數(shù)量有限，限制了公式在這些維度上的擬合能力。因此，我們假定壓降與N、M的關(guān)系是線性的。

如圖11所示，壓降與H和D的關(guān)系呈現(xiàn)非線性的高冪次關(guān)系；如圖9所示，當(dāng)其他自變量被固定，壓降隨H和D的變化趨勢都大致呈現(xiàn)二次關(guān)系。因此，我們假定壓降與H和D的關(guān)系是二次的。圖12展示了不同M下壓降隨H和D變化的散點(diǎn)圖，進(jìn)一步說明了我們假設(shè)的二次關(guān)系是合理的。

圖12 不同M下壓降隨H、D變化的散點(diǎn)圖

根據(jù)上述分析，我們提出擬合函數(shù)的形式為：

P=a×N+b×M+c2×H2+c1

×H+d2×D2+d1×D+e

(2)

式中a,b,c2,c1,d2,d1,e均為待定系數(shù)。

對參數(shù)的求解屬于回歸分析中的參數(shù)估計(jì)問題，通常采用最小二乘估計(jì)方法，最小二乘法是一種通過最小化誤差平方和尋找最佳函數(shù)匹配的優(yōu)化算法，用公式化表示為：

(3)

其中ym是擬合曲線上的點(diǎn)，yi是觀測即實(shí)驗(yàn)得到的真實(shí)值，對于一元線性回歸y=a+bx，為使上式最小可分別對a、b求偏導(dǎo)，令偏導(dǎo)等于0可計(jì)算出a、b的值?？梢酝ㄟ^數(shù)據(jù)矩陣、參數(shù)矩陣和觀測值矩陣獲得參數(shù)的顯式解，計(jì)算過程如下：

將xi,yi(i=1,2,…,n)代入y=a+bx中：

(4)

令

本文中的多元非線性回歸可以轉(zhuǎn)換成線性回歸方程求解。使用同樣的方法可以求出系數(shù)矩陣：

(5)

我們基于計(jì)算機(jī)機(jī)器學(xué)習(xí)工具庫sk-learn編程輔助計(jì)算，估計(jì)出各項(xiàng)系數(shù)的可靠經(jīng)驗(yàn)值。我們使用了472個(gè)數(shù)值模擬得到的數(shù)據(jù)并隨機(jī)分為訓(xùn)練集和測試集，其中測試集的數(shù)據(jù)占數(shù)據(jù)總數(shù)的20%。訓(xùn)練集的作用是計(jì)算擬合公式的系數(shù)，測試集的作用是對訓(xùn)練集擬合出的顯式函數(shù)結(jié)果進(jìn)行評估。通過計(jì)算，我們得到的經(jīng)驗(yàn)公式為：

(6)

式中：P——迷宮型節(jié)流件的壓降,kPa;

N——節(jié)流片的數(shù)目;

M——流過節(jié)流件的質(zhì)量流量,kg/s;

H——節(jié)流片的高度,mm;

D——相鄰節(jié)流片的間距,mm。

考慮到壓降的物理意義(不會出現(xiàn)負(fù)值)，我們對公式(3)進(jìn)行了規(guī)范，得到新的函數(shù)表達(dá)式為：

P′=Relu(P)

(7)

式中：P′——修正過的壓降;

Relu——線性整流函數(shù)，它的作用是把負(fù)數(shù)都變?yōu)?而正數(shù)保持不變。

為了對擬合出的公式精度即回歸分析的效果進(jìn)行評估，我們選擇了統(tǒng)計(jì)學(xué)中對非線性擬合常用的評估指標(biāo)決定系數(shù)R2，其含義為由自變量變化引起的因變量變化占總體因變量變化的比例，該值越大證明函數(shù)擬合性能越優(yōu)。其計(jì)算方法為：

平均觀察值:

(8)

于是可以得到總平方和:

(9)

回歸平方和和殘差平方和分別為:

(10)

(11)

最終決定系數(shù)被定義為:

(12)

我們用測試集對訓(xùn)練集擬合得到的函數(shù)進(jìn)行統(tǒng)計(jì)學(xué)評估，計(jì)算得到?jīng)Q定系數(shù)R2為0.888。在一般的統(tǒng)計(jì)學(xué)回歸分析應(yīng)用問題中，通常認(rèn)為R2大于0.8即可認(rèn)定函數(shù)擬合精度高，因此我們認(rèn)為通過回歸分析計(jì)算得到的經(jīng)驗(yàn)函數(shù)合理有效[21]。同時(shí)，擬合出的經(jīng)驗(yàn)公式通過了關(guān)于各自變量的顯著性檢驗(yàn)，說明了各自變量都對壓降有著重要的影響。此外，該經(jīng)驗(yàn)公式對數(shù)值模擬數(shù)據(jù)的均方根誤差為32.259，比對其他函數(shù)形式模型(包括N、M、H、D的高次冪模型和H、D的低次冪模型)的均方根誤差都要小，說明函數(shù)既沒有欠擬合也沒有過擬合。我們畫出了測試集真實(shí)數(shù)據(jù)和公式輸出數(shù)據(jù)的對比圖像，如圖13所示，可以較為直觀地看出擬合到的函數(shù)是可信的。

圖13 測試集真實(shí)數(shù)據(jù)和公式輸出數(shù)據(jù)的對比圖

結(jié)論：

(1)本文研究了迷宮型節(jié)流件的數(shù)值模擬方法，并通過數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對比，證明了標(biāo)準(zhǔn)k-ε模型是最適合迷宮型節(jié)流件壓降計(jì)算的湍流模型。

(2)實(shí)驗(yàn)表明，如果流通面積過小，空化現(xiàn)象將更容易出現(xiàn)，最先出現(xiàn)的位置是倒數(shù)第二塊板后側(cè)的小漩渦內(nèi)。出于安全考慮，在設(shè)計(jì)中排除能夠引起空化現(xiàn)象的幾何結(jié)構(gòu)。

(3)運(yùn)用上述經(jīng)過驗(yàn)證的湍流模型，本文研究了迷宮型節(jié)流件的壓降隨H、D、N、M的變化趨勢。當(dāng)D小于25 mm時(shí)，隨著D增加壓降顯著減??；當(dāng)D大于25 mm時(shí)，隨著D增加壓降幾乎不變。隨著H、N、M的增加，壓降都呈現(xiàn)增加的趨勢。

(4)基于大量的數(shù)值模擬數(shù)據(jù)，運(yùn)用機(jī)器學(xué)習(xí)庫sk-learn，通過多元非線性回歸方法擬合出經(jīng)驗(yàn)公式：

式中：P′——迷宮型節(jié)流件的壓降,kPa;

N——節(jié)流片的數(shù)目;

M——流過節(jié)流件的質(zhì)量流量,kg/s;

H——節(jié)流片的高度,mm;

D——相鄰節(jié)流片的間距,mm;

Relu——線性整流函數(shù)。

該公式可以被用作迷宮型節(jié)流件壓降的評估參考。

致謝

感謝中國核動力研究設(shè)計(jì)院的王嘉瑞對實(shí)驗(yàn)工作的支持與幫助。