棒束子通道两流体全隐式Picard Krylov算法

2024-01-08张宇航田兆斐李磊钱浩

哈尔滨工程大学学报 2023年12期

张宇航, 田兆斐, 李磊, 钱浩

(哈尔滨工程大学核安全与仿真技术国防重点学科实验室,黑龙江哈尔滨 150001)

两流体模型相较于均相流模型和漂移流模型有着较高的精度,而全隐式求解两流体模型能降低数值误差,以保证模型精度,但是对数值计算方法有较高的要求[1]。而目前主要的两相子通道程序,多是基于半隐式方法计算求解,如COBRAIIIC、COBRA-TF、CASTA-1、COSINE、LINDEN等两相子通道程序,存在时间步长内物理量更新不彻底引入误差的问题[2-6]。而采用全隐式算法进行子通道热工水力计算,可增强计算稳定性和计算精度。目前,学者已开发基于全隐式连续流体欧拉算法(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 两流体子通道数值计算

1.1 两流体子通道数学模型

子通道分析方法假设各流体通道利用其虚拟边界与相邻通道进行质量、能量和动量交换。该方法集总参数处理,并在控制体中利用简化横向流动的假设以完成流场求解[13-14]。子通道数学模型在传统的笛卡尔坐标系下控制方程简化,其假设为[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 两流体子通道数值离散

本文采用交错网格的方式离散子通道控制方程组,其中质量和能量方程组的控制体为标量网格,动量方程组的控制体为矢量网格。对标量网格体积积分并离散可得,全隐式液相和汽相质量离散方程为:

(9)

(10)

在式(9)和(10)中,等式左侧分别为质量随时间变化项和3个方向上流入控制体的质量流量;等式右侧分别为相变质量流量项和湍流交混和空泡漂移质量流量项。全隐式液相和汽相能量方程为:

(11)

(12)

在式(11)和(12)中,等式左侧分别为能量随时间变化项和3个方向上进入控制体的能量;等式右侧分别为压力做功项、壁面传热项、两相相变传热项、湍流交混和空泡漂移的能量项。在轴向动量网格中对轴向动量方程体积积分并离散,可得全隐式液相和汽相轴向动量离散方程为:

(13)

(14)

在式(13)和(14)中,等式左侧分别为轴向动量随时间变化项和3个方向上进入控制体的轴向动量;等式右侧分别为压力项、重力项、壁面阻力项、两相相间阻力项、轴向相变动量项。在横向动量网格中对横向动量方程体积积分并离散,可得全隐式液相和汽相横向动量离散方程(以y方向为例)为:

(15)

(16)

在式(15)和(16)中,等式左侧分别为y向动量随时间变化项、轴向进入控制体的动量、y向净流入控制体的动量;等式右侧分别为压力项、横向相变动量项、壁面阻力项、两相界面阻力项、湍流交混和空泡漂移动量项。

2 全隐式Picard Krylov方法

Picard迭代算法支持多种物理场方程隐式耦合求解,可有效减少各物理场耦合误差,且不同的物理模型可以模块化地更新,扩展性较好。而Krylov子空间算法是一种高效求解大规模稀疏线性方程组方法,使用较少的子空间将大型稀疏矩阵向量积分解,以完成高效计算。本文采用Picard迭代思想,在一个时间步长中将热工水力方程组全隐式耦合求解以保证求解精度。同时,引入Krylov子空间算法将两相流体的质量和能量方程组线性化并求解,以提升此部分的计算效率。FIPK算法流程如图1所示。在Picard迭代过程中主要包括3部分:1)热源对流场的换热计算;2)两相流体控制方程计算;3)两相流体间本构方程计算。热源对流场的换热计算需要计算棒状热源对两相流体的换热系数并进一步求解热流密度,以作为热量源相,并用于流体能量守恒方程计算。本文将动量方程和质量-能量方程分开求解,先求解动量方程以得到质量流量,用于后续的强耦合质量-能量方程的计算。最后,通过识别流型,采取对应的本构模型完成两相的界面阻力、换热面积、相变质量流量等参数计算。

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

本文将质量方程和能量方程组联立并使用Krylov子空间算法将二者耦合求解,一方面,利用Krylov算法的计算性能提升耦合方程组的求解效率;另一方面,质量和能量方程被强耦合求解,Picard迭代点列将同时满足质量和能量方程,等效于减少一组方程组,有利于提升Picard迭代收敛速度。

本文使用Krylov子空间算法中常用广义极小残量法(generalized minimum residual,GMRES)求解数值模型。由于其是一种求解线性方程组方法,所以需要对质量-能量强耦合方程组线性化,其残差方程组为:

(17)

(18)

(19)

(20)

式中:Cl、Cg分别为液相和汽相的质量残差方程组;El、Eg分别为液相和汽相的能量残差方程组。线性化质量能量耦合方程组为:

(21)

3 FIPK算法程序验证

本文使用OECD/NRC发布的《压水堆子通道及棒束通道基准题(PSBT)》中5系列含汽率分布稳态和瞬态基准题对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稳态和瞬态基准题利用X射线测量中间四通道(15、16、21、22)截面含汽率的平均值,测点位置为高部测点(3 177 mm),中部测点(2 669 mm),低部测点(2 216 mm)。

3.1 棒束通道稳态基准验证

PSBT基准题边界条件为入口温度、入口质量流量、堆芯压力和总加热功率等4个参数。为更好测试算法的适用性,本文使用边界条件范围较广的12个稳态基准题,对FIPK计算精度和计算效率验证,稳态基准题边界条件范围如表1所示。

表1 PSBT稳态棒束基准题(5系列)边界条件范围

PSBT 5稳态系列基准题数值计算结果与实验结果的空泡份额对比如图3所示,可见整体上FIPK程序在稳态计算中精度较好。仅少数算例在流速不高时计算结果在稳态含汽率计算值偏小,该误差可能与子通道横流忽略假设有关。

图3 稳态截面含汽率数值计算与实验结果对比Fig.3 Comparison chart of numerical calculation and experimental results of void fraction in steady-state section

3.2 棒束通道瞬态基准验证

本文使用PSBT 5T瞬态基准题中功率增加和流量降低2个基准题对FIPK算法瞬态验证,二者初始条件见表2[16]。

表2 PSBT瞬态棒束基准题(5T系列)初始边界条件

功率增加瞬态过程归一化边界条件如图4所示,将系统处于稳态时为计时零点,在大约第1 s时,系统近似线性提升功率,直到第5 s时功率值约为初始值的1.35倍,然后维持在此功率附近波动。该瞬态下,FIPK算法与McMaster、EDF、KAERI等机构计算结果以及PSBT实验测量值[16]对比,如图5(a)～(c)所示。由此可见,FIPK算法在功率增加瞬态计算中计算准确度较好。

图4 功率增加工况边界条件Fig.4 Boundary conditions of power increase condition

图5 功率增加工况截面含汽率随时间变化Fig.5 The void fraction changes with time under power increase condition

流量降低瞬态过程归一化边界条件如图6所示,将系统为稳态时刻为计时零点,在接近1 s时质量流量开始降低,直到约2.8 s时达到最低值约为初始值的0.55倍,然后开始回升流量直到第5 s,约为初始值的0.82倍。该瞬态下,FIPK算法与McMaster、EDF、KAERI等机构计算结果以及PSBT实验测量值[16]对比,如图7(a)～(c)所示。由此可见,FIPK在流量降低瞬态计算中计算准确度较好。