棒束子通道兩流體全隱式Picard Krylov算法

張宇航, 田兆斐, 李磊, 錢浩

(哈爾濱工程大學 核安全與仿真技術國防重點學科實驗室,黑龍江 哈爾濱 150001)

兩流體模型相較于均相流模型和漂移流模型有著較高的精度,而全隱式求解兩流體模型能降低數(shù)值誤差,以保證模型精度,但是對數(shù)值計算方法有較高的要求[1]。而目前主要的兩相子通道程序,多是基于半隱式方法計算求解,如COBRAIIIC、COBRA-TF、CASTA-1、COSINE、LINDEN等兩相子通道程序,存在時間步長內物理量更新不徹底引入誤差的問題[2-6]。而采用全隱式算法進行子通道熱工水力計算,可增強計算穩(wěn)定性和計算精度。目前,學者已開發(fā)基于全隱式連續(xù)流體歐拉算法(implicit continuous-fluid Eulerian,ICE)算法求解漂移流模型的FLICA-4、ATHAS等程序[7-9],具有較高準確度。對更精細的兩流體模型進行全隱式高效準確地計算成為了研究熱點。

目前全隱式程序廣泛采用的ICE方法本質上是一種Picard迭代算法[10]。相較于不支持全隱式的算符分解方法和求解難度較大、擴展性較差的無Jacobi矩陣的Newton Krylov算法(Jacobian-free Newton Krylov,JFNK),Picard迭代算法支持全隱式求解兩流體模型,且具有較好的擴展能力以支持模型進一步升級,但是對多物理方程組耦合問題求解效率較低[11]。而Krylov子空間算法是一種高效的大型稀疏線性方程組求解算法,能高效求解線性化的稀疏物理場方程組,能加速Picard迭代中的相關物理方程組求解,以提升Picard迭代計算效率[12]。

本文提出基于Picard迭代和Krylov子空間算法的全隱式Picard Krylov算法(fully implicit Picard Krylov, FIPK),以利用全隱式Picard算法保證計算精度,并采用Krylov子空間算法提升計算效率。

1 兩流體子通道數(shù)值計算

1.1 兩流體子通道數(shù)學模型

子通道分析方法假設各流體通道利用其虛擬邊界與相鄰通道進行質量、能量和動量交換。該方法集總參數(shù)處理,并在控制體中利用簡化橫向流動的假設以完成流場求解[13-14]。子通道數(shù)學模型在傳統(tǒng)的笛卡爾坐標系下控制方程簡化,其假設為[15]:

1) 假設流體軸向速度遠遠大于相鄰子通道間橫向速度;

2) 假設橫流從2個相鄰子通道間的虛擬間隙中流出時,此部分質量并入軸向流動方向;

3) 假設橫向的速度正交于軸向方向。

本文采用兩流體模型建模,液相質量守恒及汽相質量守恒方程分別為:

(1)

(2)

式中:αl為截面含液率;ρl為液相密度;ul、vl、wl分別為液相在x、y、z3個方向的流速;Γ為液相汽化質量流量;Tml為液相湍流交混和空泡漂移質量流量;αg為截面含汽率;ρg為汽相密度;ug、vg、wg分別為汽相在x,y,z3個方向的流速;Tmg為汽相湍流交混和空泡漂移質量流量。液相軸向動量守恒方程為:

(3)

式中:P為壓力;Fxwl為液相壁面阻力;Fxgl為汽液兩相間阻力;ΓxM為軸向液相汽化的動量。汽相軸向動量守恒方程為:

(4)

式中Fxwg為橫向液相壁面阻力。液相橫向(y方向)動量守恒方程為:

(5)

式中:Fywl為橫向液相壁面阻力;Fygl為橫向汽液兩相間阻力;ΓyM為橫向液相汽化的動量;TyMl為液相橫向湍流交混和空泡漂移的動量。汽相橫向(y方向)動量守恒方程為:

(6)

式中:Fywg為橫向汽相壁面阻力;TyMg為汽相橫向湍流交混和空泡漂移的動量。而z方向橫向動量守恒方程與y方向橫向動量守恒方程類似,本文不再贅述。液相能量守恒方程為:

(7)

式中:hl為液相比焓;qwl為壁面給液相熱流密度;Γe為兩相間換熱熱流密度;Tel為液相橫向湍流交混和空泡漂移能量。汽相能量守恒方程為:

(8)

式中:hg為汽相比焓;qwg為壁面給汽相熱流密度;Teg為汽相橫向湍流交混和空泡漂移能量。

1.2 兩流體子通道數(shù)值離散

本文采用交錯網(wǎng)格的方式離散子通道控制方程組,其中質量和能量方程組的控制體為標量網(wǎng)格,動量方程組的控制體為矢量網(wǎng)格。對標量網(wǎng)格體積積分并離散可得,全隱式液相和汽相質量離散方程為:

(9)

(10)

在式(9)和(10)中,等式左側分別為質量隨時間變化項和3個方向上流入控制體的質量流量;等式右側分別為相變質量流量項和湍流交混和空泡漂移質量流量項。全隱式液相和汽相能量方程為:

(11)

(12)

在式(11)和(12)中,等式左側分別為能量隨時間變化項和3個方向上進入控制體的能量;等式右側分別為壓力做功項、壁面?zhèn)鳠犴?、兩相相變傳熱項、湍流交混和空泡漂移的能量項。在軸向動量網(wǎng)格中對軸向動量方程體積積分并離散,可得全隱式液相和汽相軸向動量離散方程為:

(13)

(14)

在式(13)和(14)中,等式左側分別為軸向動量隨時間變化項和3個方向上進入控制體的軸向動量;等式右側分別為壓力項、重力項、壁面阻力項、兩相相間阻力項、軸向相變動量項。在橫向動量網(wǎng)格中對橫向動量方程體積積分并離散,可得全隱式液相和汽相橫向動量離散方程(以y方向為例)為:

(15)

(16)

在式(15)和(16)中,等式左側分別為y向動量隨時間變化項、軸向進入控制體的動量、y向凈流入控制體的動量;等式右側分別為壓力項、橫向相變動量項、壁面阻力項、兩相界面阻力項、湍流交混和空泡漂移動量項。

2 全隱式Picard Krylov方法

Picard迭代算法支持多種物理場方程隱式耦合求解,可有效減少各物理場耦合誤差,且不同的物理模型可以模塊化地更新,擴展性較好。而Krylov子空間算法是一種高效求解大規(guī)模稀疏線性方程組方法,使用較少的子空間將大型稀疏矩陣向量積分解,以完成高效計算。本文采用Picard迭代思想,在一個時間步長中將熱工水力方程組全隱式耦合求解以保證求解精度。同時,引入Krylov子空間算法將兩相流體的質量和能量方程組線性化并求解,以提升此部分的計算效率。FIPK算法流程如圖1所示。在Picard迭代過程中主要包括3部分:1)熱源對流場的換熱計算;2)兩相流體控制方程計算;3)兩相流體間本構方程計算。熱源對流場的換熱計算需要計算棒狀熱源對兩相流體的換熱系數(shù)并進一步求解熱流密度,以作為熱量源相,并用于流體能量守恒方程計算。本文將動量方程和質量-能量方程分開求解,先求解動量方程以得到質量流量,用于后續(xù)的強耦合質量-能量方程的計算。最后,通過識別流型,采取對應的本構模型完成兩相的界面阻力、換熱面積、相變質量流量等參數(shù)計算。

圖1 FIPK算法流程Fig.1 FIPK algorithm flow chart

本文將質量方程和能量方程組聯(lián)立并使用Krylov子空間算法將二者耦合求解,一方面,利用Krylov算法的計算性能提升耦合方程組的求解效率;另一方面,質量和能量方程被強耦合求解,Picard迭代點列將同時滿足質量和能量方程,等效于減少一組方程組,有利于提升Picard迭代收斂速度。

本文使用Krylov子空間算法中常用廣義極小殘量法(generalized minimum residual,GMRES)求解數(shù)值模型。由于其是一種求解線性方程組方法,所以需要對質量-能量強耦合方程組線性化,其殘差方程組為:

(17)

(18)

(19)

(20)

式中:Cl、Cg分別為液相和汽相的質量殘差方程組;El、Eg分別為液相和汽相的能量殘差方程組。線性化質量能量耦合方程組為:

(21)

3 FIPK算法程序驗證

本文使用OECD/NRC發(fā)布的《 壓水堆子通道及棒束通道基準題(PSBT)》中5系列含汽率分布穩(wěn)態(tài)和瞬態(tài)基準題對FIPK算法驗證。如圖2所示,B5組件加熱棒共25個,流體子通道劃分為36個,各棒源中間通道歸一化功率為1.0,邊緣通道歸一化功率為0.85。該組件軸向功率均勻分布,加熱棒長度為3 658 mm[16]。

圖2 B5組件子通道劃分示意Fig.2 B5 assembly sub-channels division diagram

PSBT 5穩(wěn)態(tài)和瞬態(tài)基準題利用X射線測量中間四通道(15、16、21、22)截面含汽率的平均值,測點位置為高部測點(3 177 mm),中部測點(2 669 mm),低部測點(2 216 mm)。

3.1 棒束通道穩(wěn)態(tài)基準驗證

PSBT基準題邊界條件為入口溫度、入口質量流量、堆芯壓力和總加熱功率等4個參數(shù)。為更好測試算法的適用性,本文使用邊界條件范圍較廣的12個穩(wěn)態(tài)基準題,對FIPK計算精度和計算效率驗證,穩(wěn)態(tài)基準題邊界條件范圍如表1所示。

表1 PSBT穩(wěn)態(tài)棒束基準題(5系列)邊界條件范圍

PSBT 5穩(wěn)態(tài)系列基準題數(shù)值計算結果與實驗結果的空泡份額對比如圖3所示,可見整體上FIPK程序在穩(wěn)態(tài)計算中精度較好。僅少數(shù)算例在流速不高時計算結果在穩(wěn)態(tài)含汽率計算值偏小,該誤差可能與子通道橫流忽略假設有關。

圖3 穩(wěn)態(tài)截面含汽率數(shù)值計算與實驗結果對比Fig.3 Comparison chart of numerical calculation and experimental results of void fraction in steady-state section

3.2 棒束通道瞬態(tài)基準驗證

本文使用PSBT 5T瞬態(tài)基準題中功率增加和流量降低2個基準題對FIPK算法瞬態(tài)驗證,二者初始條件見表2[16]。

表2 PSBT瞬態(tài)棒束基準題(5T系列)初始邊界條件

功率增加瞬態(tài)過程歸一化邊界條件如圖4所示,將系統(tǒng)處于穩(wěn)態(tài)時為計時零點,在大約第1 s時,系統(tǒng)近似線性提升功率,直到第5 s時功率值約為初始值的1.35倍,然后維持在此功率附近波動。該瞬態(tài)下,FIPK算法與McMaster、EDF、KAERI等機構計算結果以及PSBT實驗測量值[16]對比,如圖5(a)～(c)所示。由此可見,FIPK算法在功率增加瞬態(tài)計算中計算準確度較好。

圖4 功率增加工況邊界條件Fig.4 Boundary conditions of power increase condition

圖5 功率增加工況截面含汽率隨時間變化Fig.5 The void fraction changes with time under power increase condition

流量降低瞬態(tài)過程歸一化邊界條件如圖6所示,將系統(tǒng)為穩(wěn)態(tài)時刻為計時零點,在接近1 s時質量流量開始降低,直到約2.8 s時達到最低值約為初始值的0.55倍,然后開始回升流量直到第5 s,約為初始值的0.82倍。該瞬態(tài)下,FIPK算法與McMaster、EDF、KAERI等機構計算結果以及PSBT實驗測量值[16]對比,如圖7(a)～(c)所示。由此可見,FIPK在流量降低瞬態(tài)計算中計算準確度較好。

圖6 流量降低瞬態(tài)工況各時刻邊界條件Fig.6 Boundary conditions of flow reduction condition

圖7 流量降低工況截面含汽率隨時間的變化Fig.7 The void fraction changes with time under flow reduction condition

3.3 穩(wěn)態(tài)和瞬態(tài)工況計算效率

FIPK算法使用Krylov子空間算法將質量和能量方程強耦合求解,能夠有效地提升計算效率。本文分別使用FIPK算法和傳統(tǒng)Picard算法對PSBT穩(wěn)態(tài)和瞬態(tài)基準題求解,并記錄計算用時。將計算用時的倒數(shù)作為計算效率,并將FIPK算法與Picard算法計算效率之商作為二者的相對計算效率。

在硬件配置為32 GB內存、Inter i7-8700 CPU環(huán)境下,穩(wěn)態(tài)PSBT基準題求解中,FIPK相對Picard計算效率如表3所示。FIPK相較于Picard,在穩(wěn)態(tài)基準題計算中計算效率最高提升89.74%;在功率增加和流量降低的瞬態(tài)基準題計算中,計算效率分別提升1.29%和20.68%。

表3 FIPK與Picard算法計算效率對比

其中,少數(shù)基準題計算中,FIPK算法相較于傳統(tǒng)Picard算法計算效率提升不明顯,可能因為:1)在前后時刻迭代解相差不大時,數(shù)值計算過程收斂較容易,FIPK算法與傳統(tǒng)Picard算法計算效率相當;2)少數(shù)情況下,質量和能量方程的求解效率可能不是影響整體計算效率的主要因素。

4 結論

1)本文基于Picard算法和Krylov子空間算法,提出的一種全隱式兩流體子通道計算方法FIPK,能準確計算兩相棒束熱工水力穩(wěn)態(tài)和瞬態(tài)現(xiàn)象。

2)FIPK算法采用全隱式Picard算法實現(xiàn)導熱方程、對流擴散方程、兩相本構方程緊耦合求解,保證計算精度。同時,采用Krylov子空間算法實現(xiàn)質量和能量方程強耦合求解,提升了計算效率。

3)FIPK算法通過了PSBT穩(wěn)態(tài)基準題、功率提升和流量降低瞬態(tài)基準題的驗證,數(shù)值結果表明FIPK算法具有較高的精度和效率。