基于JFNK的兩流體三流場兩相流模型的全隱算法研究

樊 杰 吳和鑫 茍軍利 單建強(qiáng)

（西安交通大學(xué)核科學(xué)與技術(shù)學(xué)院 西安 710049）

對(duì)于環(huán)狀流下干涸以及干涸后的傳熱以及核反應(yīng)堆再淹沒傳熱，液滴的傳熱傳質(zhì)在其中扮演了重要的角色，目前主流的核反應(yīng)堆系統(tǒng)分析程序大都采用兩流體兩相流模型，沒有模擬液滴效應(yīng)，傳熱計(jì)算有一定的誤差，而三流場兩相流模型考慮了液滴的流動(dòng)傳熱，成為核反應(yīng)堆安全分析程序發(fā)展的方向。

目前國際上基于兩流體三流場的兩相流模型或者 核 反 應(yīng) 堆 安 全 分 析 程 序 有CATHARE3［1］、SPACE［2］、CUPID［3］、TRAC-M［4］、C3DM［5］以 及Frepoli等［6］和筆者所在的核安全與運(yùn)行實(shí)驗(yàn)室所開發(fā)的三流場兩相流模型［7］，除Frepoli開發(fā)的基于漂移流的兩流體三流場模型外，其余這些模型或程序都采用半隱數(shù)值算法，與其他基于半隱數(shù)值算法的核反應(yīng)堆安全分析程序一樣，例如：RELAP5，數(shù)值穩(wěn)定性受聲速庫侖特值的影響。而全隱數(shù)值算法是絕對(duì)穩(wěn)定的，因此也成為核反應(yīng)堆系統(tǒng)分析研究的一個(gè)熱點(diǎn)。

目前大部分的基于全隱數(shù)值算法的程序都基于傳統(tǒng)的牛頓迭代法，例如：Frepoli開發(fā)的兩流體三流場模型、CATHARE2［8］等，而傳統(tǒng)的牛頓迭代法需要形成雅克比矩陣，對(duì)于本構(gòu)模型較多的兩流體模型雅克比矩陣的形成具有一定的困難，因此近年來一種無需雅克比矩陣的牛頓-Krylov迭代法（Jacobianfree Newton-Krylo，JFNK）被用于兩相流模型的全隱數(shù)值算法，當(dāng)前JFNK用于兩相流模型全隱數(shù)值算法的研究還不夠成熟，且集中于漂移流模型或者兩流體模型，如：Ashrafizadeh等［9］、Hajizadeh等［10］、Hu等［11］、Wu等［12］的漂移流模型，Mousseau等［13］、Zou等［14］的兩流體模型，在公開文獻(xiàn)中基于JFNK算法的兩流體三流場模型的全隱算法的相關(guān)研究還沒有，因此本研究嘗試對(duì)先進(jìn)的兩流體三流場兩相流模型開發(fā)基于JFNK算法的全隱數(shù)值算法。

本研究針對(duì)兩流體三流場兩相流模型，建立基于JFNK算法的全隱算法，并針對(duì)三流場中某一相或兩相缺失造成的數(shù)值計(jì)算困難進(jìn)行了特殊的數(shù)值處理。程序模擬了Ransom水龍頭數(shù)值標(biāo)準(zhǔn)題，以驗(yàn)證數(shù)值算法的精度，并通過模擬過冷沸騰實(shí)驗(yàn)和Dryout傳熱實(shí)驗(yàn)以驗(yàn)證程序計(jì)算單相、單相到兩相的過渡、沒有液滴的兩相和三流場兩相等工況下數(shù)值算法的精度。

1 數(shù)學(xué)模型

兩相流數(shù)學(xué)模型采用一維兩流體三流場模型，分別對(duì)氣相、連續(xù)液相和彌散液滴相建立質(zhì)量、動(dòng)量和能量方程，模擬了氣液之間的傳熱傳質(zhì)以及連續(xù)液相與彌散液滴相之間的夾帶沉降，氣相的質(zhì)量守恒方程為：

式中：右邊項(xiàng)分別為單位體積的相間傳質(zhì)量和壁面?zhèn)髻|(zhì)量。

液相的質(zhì)量守恒方程為：

式中：右邊項(xiàng)分別為連續(xù)液相與氣相的相間換熱量、壁面?zhèn)髻|(zhì)量、液滴夾帶量、液滴沉降量。

液滴相的質(zhì)量守恒方程為：

氣相的動(dòng)量守恒方程為：

式中：右邊項(xiàng)分別為彌散液滴與氣相的相間換熱量、壁面?zhèn)髻|(zhì)量、液滴夾帶量、液滴沉降量。

式中：右邊項(xiàng)分別為壓力梯度項(xiàng)、重力項(xiàng)、壁面摩擦力、局部阻力、相間傳質(zhì)引起的動(dòng)量轉(zhuǎn)移、氣相與連續(xù)液相的相間摩擦、氣相與彌散液滴的相間摩擦。

液相的動(dòng)量守恒方程為：

式中：右邊項(xiàng)分別為壓力梯度項(xiàng)、重力項(xiàng)、壁面摩擦力、局部阻力、氣相與連續(xù)液相的相間摩擦、相間傳質(zhì)引起的動(dòng)量轉(zhuǎn)移、液滴夾帶沉降引起的連續(xù)液相的動(dòng)量轉(zhuǎn)移。

液滴相的動(dòng)量守恒方程為：

式中：右邊項(xiàng)分別為壓力梯度項(xiàng)、重力項(xiàng)、相間傳質(zhì)引起的動(dòng)量轉(zhuǎn)移、液滴夾帶沉降引起的彌散液滴的動(dòng)量轉(zhuǎn)移、氣相與彌散液滴的相間摩擦。

氣相的能量守恒方程為：

連續(xù)液相和彌散液滴相被認(rèn)為是熱平衡的，液相和液滴相合成一個(gè)液相的能量方程，即：

式中：右邊項(xiàng)分別為壓力做功項(xiàng)、壁面換熱量、相間換熱量、傳質(zhì)引起的能量轉(zhuǎn)移項(xiàng)。

兩流體模型守恒方程封閉需要的本構(gòu)模型包括相間換熱、相間摩擦、壁面換熱、壁面摩擦、液滴夾帶沉降等模型，這些本構(gòu)模型與流型和換熱模式有關(guān)，程序中根據(jù)計(jì)算工況的需要選用相應(yīng)流型或換熱模式下的模型，各模型的具體形式詳見文獻(xiàn)［7］。

2 數(shù)值算法

2.1 全隱數(shù)值離散

守恒方程離散采用基于交錯(cuò)網(wǎng)格的有限體積法，時(shí)間差分采用向后的歐拉法，空間離散采用Van Albada［15］高階精度差分格式，求解變量依次為氣相內(nèi)能、液相內(nèi)能、空泡份額、液相份額、壓力、氣相速度、液相速度和液滴速度。

對(duì)于控制體i，其守恒方程有汽相能量方程FUg，i、液相能量方程FUl，i、汽相質(zhì)量方程Fαg，i、連續(xù)液相質(zhì)量方程Fαf，i、彌散液滴質(zhì)量方程FPi，對(duì)于接管i+1/2，其守恒方程有汽相動(dòng)量方程FVg，i+12、液相動(dòng)量方程FVf，i+12、液滴相動(dòng)量方程FVd，i+12。此處展示了彌散液滴相質(zhì)量守恒方程、汽相動(dòng)量守恒方程、液相能量守恒方程等3個(gè)方程的全隱式離散形式。

彌散液滴相質(zhì)量守恒方程：

汽相動(dòng)量守恒方程：

液相能量守恒方程：

2.2 JFNK算法

對(duì)于由所有控制體質(zhì)量和能量離散方程以及所有接管動(dòng)量離散方程形成的方程組，由JFNK算法求解，外迭代是經(jīng)典牛頓法，牛頓迭代的基本格式為：

對(duì)于該線性方程組，使用Krylov子空間迭代法進(jìn)行求解，常用的Krylov子空間迭代法為GMRES算法。迭代過程中無需寫出雅克比矩陣，只需要雅克比矩陣與向量的點(diǎn)積，可以使用差分近似代替微分，即：

式中：ε為差分步長或者擾動(dòng)參數(shù)［16］：

式中：N為線性方程組中方程的個(gè)數(shù)；a為機(jī)器精度。

對(duì)于JFNK算法，合適的預(yù)處理技術(shù)可以有效地加快Krylov子空間迭代的收斂。JFNK方法中常用的線性預(yù)處理方式有：左預(yù)處理、中間預(yù)處理、右預(yù)處理等，由于右預(yù)處理技術(shù)不會(huì)改變方程的右邊殘差項(xiàng)，有利于JFNK方法中Krylov子空間迭代收斂準(zhǔn)則的選取，因此本研究中采用右預(yù)處理，即：

預(yù)處理矩陣P與雅克比矩陣J越接近，Krylov子空間迭代收斂越快，由于半隱式離散得到的雅克比矩陣與全隱式離散得到的雅克比矩陣在形式和量級(jí)方面都是相似的，因此半隱式離散的雅克比矩陣被當(dāng)作預(yù)處理矩陣。

2.3 單相問題數(shù)值處理

當(dāng)使用兩流體三流場模型計(jì)算兩相流時(shí)，可能會(huì)遇到某一相或者兩相不存在，例如：單相水、液滴不存在的泡狀流和彈狀流、連續(xù)液相不存在的彌散流、單相汽等，缺失相的守恒方程是不存在的，此時(shí)線性方程組和預(yù)處理矩陣都是奇異的，之前基于半隱算法三流場模型的研究中，給缺失相一個(gè)很小的份額用于計(jì)算相間換熱和相間摩擦，使得相間換熱和相間摩擦不為零，可以解決線性方程組奇異的問題。例如：對(duì)于泡狀流，由于沒有液滴，因此液滴的動(dòng)量守恒方程簡化為：

由于使用很小的液滴份額計(jì)算相界面面積，因此相間摩擦系數(shù)不為零，這樣得到的液滴速度等于汽相速度。然而，缺失相的份額在守恒方程中仍然是0，當(dāng)進(jìn)行“Jacobian-Free”步驟時(shí)，新的(x k+ε·x k)中液滴份額仍然是0，在差分步中液滴份額沒有貢獻(xiàn)，因此會(huì)產(chǎn)生非物理的解。

因此在當(dāng)前研究中，汽相和連續(xù)液滴相的份額為10-5~0.999 99，液滴份額為10-4~0.999 9，經(jīng)過這樣的數(shù)值處理，相間摩擦和相間換熱是連續(xù)的，在各種流型之間轉(zhuǎn)換時(shí)不存在數(shù)值困難。且在“Jacobian-Free”步驟中，(xk+ε·x k)中各相的份額均不等于x k中的相份額，因此不會(huì)產(chǎn)生非物理的解。雖然缺失相的源項(xiàng)不等于0，但相份額很小，數(shù)值誤差時(shí)是可以接受的。

2.4 程序計(jì)算流程

程序的計(jì)算流程如圖1所示，主要包括JFNK算法模塊和本構(gòu)模型計(jì)算模塊。JFNK算法模塊包括矩陣預(yù)處理、GMRES迭代和牛頓迭代，使用JFNK完成線性方程組的求解后進(jìn)入下一時(shí)間步的計(jì)算直至瞬態(tài)計(jì)算結(jié)束。在JFNK算法模塊中，x以及離散方程組的殘差F（x）在不斷地更新，需要通過本構(gòu)模型計(jì)算模塊不斷地賦值計(jì)算。本構(gòu)模型計(jì)算模塊首先需要通過x k中的壓力和內(nèi)能計(jì)算其他物性參數(shù)，然后計(jì)算相間換熱、相間摩擦、壁面換熱等守恒方程中其他項(xiàng)。

圖1 程序計(jì)算流程Fig.1 Flowchartofthecodecalculation

3 程序驗(yàn)證

3.1 水龍頭問題

為了驗(yàn)證數(shù)值算法的穩(wěn)定和精度，數(shù)值實(shí)驗(yàn)基準(zhǔn)題——Ransom水龍頭問題［17］被模擬，實(shí)驗(yàn)的初始和邊界條件如表1所示。

表1 實(shí)驗(yàn)初始條件和邊界條件Table1 Initialandboundaryconditionofthewaterfaucet test

研究對(duì)象為12 m長的豎直向上的圓管，直徑1 m，初始時(shí)刻圓管內(nèi)氣相份額為0.2，氣相速度為0.0 m?s-1，液相速度為10 m?s-1，液相溫度為50℃，壓力為105Pa。在瞬態(tài)過程中，進(jìn)口的空泡分?jǐn)?shù)、液相和氣相速度保持不變，出口壓力恒定105Pa不變。瞬態(tài)過程如圖2所示，初始液膜與管壁平行，隨著瞬態(tài)的進(jìn)行，由于重力的作用，液相加速下降。在受進(jìn)口邊界條件影響的區(qū)域，沿著管壁液相份額不斷地減少，而不受進(jìn)口條件影響的區(qū)域，液相做自由落體運(yùn)動(dòng)，空泡份額一直保持不變，因此沿著管壁空泡份額不連續(xù)。最終液相和氣相形成穩(wěn)態(tài)。

圖2 水龍頭問題示意圖Fig.2 Sketchofthewaterfaucettest

空泡份額和液相速度的解析解為：

圖3展示了瞬態(tài)開始后0.5s時(shí)程序計(jì)算的空泡份額和液相速度與解析解和NUSOL-SYS程序（西安交通大學(xué)核安全與運(yùn)行實(shí)驗(yàn)室開發(fā)的系統(tǒng)分析程序）計(jì)算結(jié)果的比較?？刂企w長度為0.3m，時(shí)間步長為0.005s，程序計(jì)算值與實(shí)驗(yàn)值符合得很好。由于空間離散使用了VanAlbada高階精度的差分格式，較NUSOL-SYS使用的一階迎風(fēng)差分格式，數(shù)值耗散更小，因此程序的計(jì)算值比NUSOL-SYS程序的計(jì)算值更加精確。比較結(jié)果證明了當(dāng)前全隱數(shù)值算法的精度比較高，且由于0.5s時(shí)空泡份額的分布在0.2~0.5，此時(shí)沒有液滴夾帶產(chǎn)生，該工況的精確模擬也論證了對(duì)于某一相或兩相缺失的問題而進(jìn)行的數(shù)值處理是合適的。

圖3 0.5 s時(shí)全隱算法計(jì)算結(jié)果與解析解和RELAP5結(jié)果的比較 (a)空泡份額，(b)液相速度Fig.3 Comparison between the predicted results by fully-implicit numerical algorithm,analytical solution and the predicted results by semi-implicit numerical algorithm at 0.5 s (a)Void fraction,(b)Liquid velocity

3.2 Bartolomei過冷沸騰實(shí)驗(yàn)

為了評(píng)估程序計(jì)算單相傳熱以及單相到兩相過渡時(shí)的表現(xiàn)，這部分模擬了Bartolomei單管過冷沸騰傳熱實(shí)驗(yàn)［18］，模擬選取了中壓中流速中過冷的工況2、中壓低流速低過冷的工況5、中壓高流速高過冷的工況8以及一組低壓的工況22，實(shí)驗(yàn)條件如表2所示，為了平衡總的計(jì)算時(shí)間和Krylov子空間迭代的收斂速度，使得總的計(jì)算時(shí)間盡可能少，各工況的CFL（CFL=最大相速度×?xí)r間步長/網(wǎng)格長度，該值越大，Krylov子空間迭代收斂越慢）值為0.276~0.455。

表2 Bartolomei過冷沸騰實(shí)驗(yàn)條件Table 2 Experimental conditions of Bartolomei subcooledboiling experiment

當(dāng)具有一定過冷度的水進(jìn)入豎直向上的圓管后，沿著流動(dòng)方向不斷地被加熱，當(dāng)達(dá)到過冷沸騰起始點(diǎn)后過冷沸騰發(fā)生，沿著流動(dòng)方向空泡份額的分布如圖4所示?；谌[算法的預(yù)測結(jié)果與實(shí)驗(yàn)值符合的較好，且與之前開發(fā)的基于半隱算法的兩流體三流場模型的計(jì)算結(jié)果基本一致，同時(shí)由于圖4中計(jì)算的幾組過冷沸騰工況均沒有達(dá)到環(huán)霧流，即沒有液滴，而兩流體三流場模型在單相、泡狀流和彈狀流的模型與NUSOL-SYS程序的模型一致，故基于全隱算法的模擬結(jié)果也與NUSOL-SYS程序的計(jì)算結(jié)果符合的很好。此外，結(jié)果也論證了當(dāng)前模型和數(shù)值算法在單相、從單相到兩相的過渡以及缺少液滴的兩相工況都表現(xiàn)的很好。

圖4 全隱算法計(jì)算結(jié)果與實(shí)驗(yàn)值和NUSOL-SYS、半隱算法結(jié)果的比較 (a)工況2，(b)工況5，(c)工況8，(d)工況22Fig.4 The comparison between the predicted results by fully-implicit numerical algorithm,experimental results and the predicted results by NUSOL-SYS and semi-implicit numerical algorithm at 0.5 s (a)Test 2,(b)Test 5,(c)Test 8,(d)Test 22

3.3 Dryout傳熱實(shí)驗(yàn)

為了驗(yàn)證基于全隱數(shù)值算法的兩流體三流場兩相流模型計(jì)算彈狀流、環(huán)霧流、彌散流之間過渡時(shí)的表現(xiàn)，Dryout傳熱實(shí)驗(yàn)被模擬。在環(huán)狀流區(qū)域，由于氣芯不斷地剪切液膜而夾帶產(chǎn)生液滴，液膜厚度不斷地變小，當(dāng)液膜厚度減小到臨界值以下時(shí)，流體不能帶走足夠的熱量，壁面溫度飛升，發(fā)生沸騰臨界現(xiàn)象。模擬選取了Bennett［19］和Becker［20］的實(shí)驗(yàn)中兩種比較極端的工況，即干涸發(fā)生時(shí)液滴較多的高流量工況Becker-221、Becker-232和干涸發(fā)生時(shí)液滴較少的低流量工況Becker-277、Bennett-5358，實(shí)驗(yàn)條件見表3。

表3 Dryout傳熱實(shí)驗(yàn)條件Table 3 Experimental conditions of dryout experiment

圖5展示了4組實(shí)驗(yàn)工況沿軸向汽相份額和液滴份額的分布。對(duì)于工況Becker-277和Bennett-5358，干涸點(diǎn)對(duì)應(yīng)的汽相份額大于0.99，干涸后區(qū)域液滴份額較少，液滴對(duì)換熱的影響很小，干涸后壁面與蒸汽的對(duì)流換熱占據(jù)壁面?zhèn)鳠嶂饕糠帧?/p>

圖5 汽相份額和液滴份額的軸向分布Fig.5 The axial fraction distribution of void and droplet

圖6顯示了基于半隱數(shù)值算法和全隱數(shù)值算法的三流場模型計(jì)算得到的壁溫與NUSOL-SYS程序計(jì)算得到的壁溫基本相同，且全隱數(shù)值算法的計(jì)算結(jié)果與半隱數(shù)值算法的計(jì)算結(jié)果基本一致。而對(duì)于高流量的工況（Becker-221和Bennett-232），干涸發(fā)生時(shí)汽相份額約為0.9，液滴份額較多，由于彌散的液滴蒸發(fā)吸收了很多熱量，壁面溫度不會(huì)一直上升。圖7中計(jì)算值與實(shí)驗(yàn)值比較可以發(fā)現(xiàn)，由于液滴的作用，三流場模型計(jì)算的壁溫比NUSOL-SYS程序計(jì)算的壁溫更加準(zhǔn)確，且液滴的夾帶沉降模型以及液滴的相間換熱模型是精確的，同樣全隱數(shù)值算法的計(jì)算結(jié)果與半隱算法的計(jì)算結(jié)果基本相同。

圖6 低流量全隱數(shù)值算法的計(jì)算結(jié)果與實(shí)驗(yàn)值、NUSOLSYS計(jì)算值和半隱數(shù)值算法的計(jì)算結(jié)果比較Fig.6 The comparison between the predicted results by fullyimplicit numerical algorithm,experimental results,simulated results by NUSOL-SYS and the predicted results by semiimplicit numerical algorithm for the conditions with low mass flux

圖7 高流量全隱數(shù)值算法的計(jì)算結(jié)果與實(shí)驗(yàn)值、NUSOLSYS計(jì)算值和半隱數(shù)值算法的計(jì)算結(jié)果比較Fig.7 Comparison between the predicted results by fullyimplicit numerical algorithm,experimental results,simulated results by NUSOL-SYS and the predicted results by semiimplicit numerical algorithm for the conditions with high mass flux

干涸后傳熱實(shí)驗(yàn)的模擬結(jié)果也說明了當(dāng)前開發(fā)的兩流體三流場模型的全隱數(shù)值算法在計(jì)算彈狀流-環(huán)霧流、環(huán)霧流-彌散流之間過渡時(shí)表現(xiàn)很好。

4 結(jié)語

本文針對(duì)兩流體三流場兩相流模型，實(shí)現(xiàn)了基于JFNK算法的全隱數(shù)值算法，并通過模擬相關(guān)實(shí)驗(yàn)對(duì)數(shù)值算法進(jìn)行驗(yàn)證，得出以下結(jié)論：

1）Ransom水龍頭數(shù)值基準(zhǔn)題的模擬結(jié)果表明：基于Van Albada高階精度差分格式的全隱數(shù)值算法比NUSOL-SYS程序具有更高的精度，且當(dāng)前數(shù)值算法對(duì)單相、沒有液滴相的兩相數(shù)值處理是可行的。

2）過冷沸騰實(shí)驗(yàn)的模擬結(jié)果表明程序可以成功從單相水過渡到彈狀流，空泡份額分布與實(shí)驗(yàn)值符合較好，且由于過冷沸騰實(shí)驗(yàn)中沒有液滴存在，程序計(jì)算結(jié)果準(zhǔn)確性與NUSOL-SYS程序的計(jì)算結(jié)果相當(dāng)。同時(shí)當(dāng)前全隱數(shù)值算法在單相到兩相過渡以及泡狀流-彈狀流過渡時(shí)表現(xiàn)很好。

3）Dryout傳熱實(shí)驗(yàn)的模擬結(jié)果表明，兩流體三流場模型計(jì)算的壁溫比傳統(tǒng)的兩流體模型更加準(zhǔn)確，且當(dāng)前全隱算法的計(jì)算結(jié)果接近半隱數(shù)值算法的精度。數(shù)值算法計(jì)算彈狀流-環(huán)霧流、環(huán)霧流-彌散流之間過渡時(shí)表現(xiàn)很好。

本研究開發(fā)的全隱數(shù)值算法的核反應(yīng)堆安全分析程序的發(fā)展奠定基礎(chǔ)，但為了提高系統(tǒng)分析程序的模擬精度，還需要進(jìn)一步的研究來改善兩相流中一些模型的計(jì)算精度，如CHF模型、一些相間傳熱模型。

作者貢獻(xiàn)聲明樊杰：模型和算法開發(fā)，程序設(shè)計(jì)，程序驗(yàn)證，文章撰寫；吳和鑫：實(shí)驗(yàn)數(shù)據(jù)調(diào)研和處理，文章修改；茍軍利：論文整體設(shè)計(jì)，對(duì)文章作批評(píng)性審閱，研究經(jīng)費(fèi)支持；單建強(qiáng)：提供技術(shù)指導(dǎo)，對(duì)文章作批評(píng)性審閱。