高温气冷堆换料过程的球流数值模拟与分析

2022-09-11严安孙喜明董玉杰

哈尔滨工程大学学报 2022年8期

严安，孙喜明，董玉杰

(清华大学核能与新能源技术研究院，北京 100084)

高温气冷堆具有发电效率高、安全性高、用途广泛等优点，是符合第4代先进核能系统技术要求的堆型之一[1]。我国正在建设的250 MW模块式高温气冷堆(HTR-PM)采用陶瓷型包覆颗粒弥散在石墨基体中制成的球形燃料元件，球床堆芯由统一直径的球形燃料元件在重力作用下随机自然堆积形成。球床堆在运行过程中不停堆换料，燃料元件在堆芯中多次循环达到一定燃耗深度后卸出堆芯。堆芯中的球流运动是一种典型的稠密颗粒流，由于球床运动呈现间歇性特点，球床大多数时间处于静力平衡状态，球体间存在长时间力链网络，惯性可以忽略不计，又属于准静态流动。球床堆芯的行为对堆芯中子物理计算、热工水力性能和安全分析有直接影响[2-3]，有必要对这种准静态稠密颗粒流的关键特性进行研究。

稠密颗粒流广泛存在于化学工程、岩土工程等领域，目前尚不存在一套完整的连续性方法可用于描述不同几何边界下的平均流动特性，因此离散单元方法(discrete element method, DEM)由于可以获得颗粒流动的微观细节而得到广泛应用[4]。Li等[5]应用DEM对二维料斗中的双粒径颗粒流进行研究，分析了局部堆积结构对整体流率的影响。Masson等[6]使用DEM分析了料斗中壁面压力的分布。DEM模拟被逐步引入至球床堆领域的研究中。Yang等[7-8]模拟了10兆瓦高温气冷堆(HTR-10)三维球床堆芯的静态堆积特性，重点分析了孔隙率在堆芯中的分布规律。Rycroft等[9]通过对44万个球体的球床模拟研究，考察了球流平均速度、搅浑和扩散效应以及滞留时间的分布规律。Li等[10]通过考察摩擦系数对于HTR-10堆芯流型和滞留区的影响，认为球体-壁面摩擦系数比球体-球体摩擦系数起到更显著的作用。Jia等[11]针对假定的60°锥角球床堆芯切片进行模拟，分析了改变换料速率对于球床流速分布和流型的影响。Gui等[12]研究了HTR-PM堆芯几何边界的变化对球体滞留与流线分布规律的影响。Qi等[13]模拟了HTR-10堆芯切片的堆积和换料，得出了弹性恢复系数和球流均匀性分布之间的关系。通常认为，DEM模拟中球体刚度越大需要的计算量越大。高温堆堆芯中燃料球石墨基体的杨氏模量在1010Pa左右，为了提高计算效率，上述研究均不同程度减小了球体杨氏模量的取值，一般在106～5×108Pa。Sun等[14]针对德国HTR-MODUL球床堆芯的换料过程的准静态球流进行了DEM模拟，由于采用了接近真实的球体刚度，模拟耗费机时超过半年。对于准静态颗粒流动过程，球体刚度是影响流动特性的主要参数[15]，这对于球床式高温堆的换料过程至关重要。上述文献没有就刚度对换料球流的影响展开研究，因此有必要在以往研究的基础上，进一步探究不同球体刚度取值下的球床模拟结果与真实形态的相关程度。

本文应用DEM，模拟了高温气冷堆球床堆芯的准静态流动过程，分析了杨氏模量参数对于换料过程的影响，研究了保证球流DEM计算不失真的最低刚度值；分析了不同球体数量的球床在循环换料过程中的球流运动特性以及球床在换料过程中流型转变与滞留现象。

1 建模方法与球床堆芯结构

1.1 建模方法

DEM是一种模拟非连续性物质的数值方法，一般分为软球模型和硬球模型[3]。鉴于硬球模型无法实现多体接触模拟，本文基于软球模型原理，将堆芯中球型燃料元件视为相互独立的非完全刚性球体单元，相互致密接触并作用，且因接触作用使球体间接触处产生形变，导致球体相互间发生区域变形。换料准静态运动中各球体呈现空间运动，根据动量与动量矩定理，第i个球体的基本运动方程为：

(1)

(2)

球体之间的接触力可分为法向力和切向力：

Fc=Fnn+Ftt

(3)

Fn=f(δ/d)(knδ-γnUn)

(4)

Ft=f(δ/d)(ktSt-γtUt)

(5)

式中：d为球体直径；kn、kt分别为法向和切向的刚度；δ、St分别为球体间重叠量和弹性切向位移，弹性切向位移在接触发生时刻为0；γn、γt分别为法向和切向的阻尼系数；Un、Ut分别为法向和切向的相对速度分量。

(6)

式中：E*为两接触球体的等效杨氏模量；R*为两接触球体的等效半径。

DEM方法建立的数值模型通常采用中心差分、显示时间积分格式求解，是典型的刚性系统，收敛条件由时间步长决定。根据文献[17]，对于单自由度系统，即一个球体单元而言，临界时间步长为：

(7)

可见DEM模拟中使用的杨氏模量越大，则临界时间步长越小，完成相同模拟所需时间步越多，计算效率越低。

1.2 球床堆芯结构

本文的研究对象为HTR-PM的球床堆芯，HTR-PM的一、二回路结构如图1所示。反应堆堆芯中大约有4.2×105个燃料球，燃料球直径d为0.06 m，形成了直径3 m、高11 m的堆芯，堆芯周边的反射层是耐高温的石墨[1]。堆芯底部锥角为30°，卸料管直径为0.5 m。图2是正常运行工况下换料过程的堆芯底锥部分结构示意图。燃料球在堆芯中向下运动直至进入卸料管，堆芯顶部加入等数目的再循环燃料球或者新燃料球。

图1 模块式高温气冷堆结构示意Fig.1 Structure schematic diagram of the HTR-PM

图2 球床堆芯底部锥体结构示意Fig.2 Schematic diagram of cone area in pebble bed

1.3 模型参数选取及换料过程模拟

文献[14]针对德国HTR-MODUL球床堆芯换料过程进行了DEM模拟，其换料过程、球体单元和堆芯几何等设计参数与HTR-PM类似，文献中选取的材料物性参数已经过验证，因此本文的物性参数与其保持一致，相应DEM模拟所使用的主要参数如表1所示。其中，密度为燃料球的真实属性值；泊松比、恢复系数和摩擦系数直接影响接触力的计算。本文模拟中共采用3种DEM模拟常见的杨氏模量，分别为9.8×105Pa、9.8×106Pa、9.8×107Pa，探究杨氏模量对于换料过程球流的影响。本文分别使用8万、10万、15万、20万个燃料球，堆积形成高度为2.3、2.8、4.0、5.1 m的球床堆芯(高度是以卸料口为起始位置，至球床堆芯顶部)，模拟堆芯高度对于换料过程的影响。上述球体总数也包含卸料管中的球体。需要指出的是，堆芯及卸料管壁面的杨氏模量固定取为9.8×109Pa。球-球与球-壁面之间的接触参数相同，参见表1。

表1 模型主要参数Table 1 Main parameters used in model

换料过程模拟分为2步：1)球体在重力作用下堆积，直至球床高度不再变化，即堆芯进入稳定状态；2)卸料管末端区域按固定的频率删除球体单元，同时，堆芯顶部填入相同数量的球体单元。真实工程下，HTR-PM换料速率约为4 个/min。为加速DEM计算，已知文献选取的换料速率均远高于该值，一般为50～200 个/s。为更准确地模拟间歇性换料过程，本文模拟中球体的换料速率为8 个/s，为已知文献中最接近真实工况的换料速度。

2 换料模拟结果与分析

2.1 刚度对于球床换料过程的影响

图3为不同刚度模拟中3个时刻的球床切片，其中图3(a)～(c)为刚度取9.8×105Pa的球床，图3(d)～(f)为刚度取9.8×106Pa的球床，图3(g)～(i)为刚度取9.8×107Pa的球床；图3(a)、(d)、(g)为初始时刻t=0 s，图3(b)、(e)、(h)为中间时刻t=500 s，图3(c)、(f)、(i)为最终时刻t=1 000 s。对比图3(a)、(d)、 (g)可知，由于模拟中刚度的取值相差100倍，球床的整体高度明显不同。为直观显示球床运动的差异，在初始时刻，分别将3种球床每6d高度处的2层球体标记为不同颜色。可以看出，随着换料过程的持续，球床逐渐向下运动，球流的分布出现明显差异，图3(f)、(i)与图3(c)相比，靠近壁面的区域流动较为缓慢，呈现出更突出的深“V”型。

图3 不同杨氏模量换料过程球床球流切面Fig.3 Cross section snapshots of pebble flow for pebble beds with different Young′s modulus

为了进一步定量分析球流的运动规律，需要分析球床的速度场。由于球体间歇性卸出，仅小部分球体会发生运动，球床大多数时间为静止状态，很难获得球体的速度信息，所以本文利用同一球体在相邻卸料时间步(即t0和t0+Δt)时的位置差值与时间的比值，作为球体的平均运动速度予以分析。

图4为换料过程中平均轴向(z向)速度分布，分别取球床高度H=15d，20d，24d处进行定量分析。可以看出，由于卸料管和底锥的影响，堆芯中心球体的轴向流速相比壁面附近的更高，且随着球体越贴近卸料口，即球体高度越低，这一现象越明显。低刚度(E=9.8×105Pa)的模拟中，球体轴向速度在径向上的分布整体相对平均，同时不同高度的流速差别较小；但随着刚度的增加，球床下部(z=15d、20d)半径大于20d处靠近壁面的球体轴向速度趋近于零，发生了明显的滞留，球床上部的球流速度相对平缓。

图4 不同杨氏模量与高度下球流轴向速度分布Fig.4 Velocity profiles for pebble flow at several cross sections for beds with different Young′s modulus

图5 不同杨氏模量下球床球流流线Fig.5 Streamline profiles for pebble flow with different Young′s modulus

不同球体刚度换料过程模拟的球床流线如图5所示，球流速度矢量如图6所示。选取速度统计中的位移中点作为每个速度矢量的起点，速度矢量的方向指向球体位移的方向，大小和颜色由速度大小指定。不难看出，流线在卸料口汇聚，球流在底锥区域内发生搅混，这说明球体受重力的驱动，在壁面支撑力、摩擦力与其他球体接触力的共同作用下，互相竞争排出堆芯。

图6 不同杨氏模量下球床球流速度矢量Fig.6 Velocity vector profiles for pebble flow with different Young′s modulus

由图5(a)可以看出，低刚度的球床在换料过程中不存在滞留，流动状态趋向于整体流，高度大于15d的部分没有明显的径向运动。由图5(b)、(c)可以看出，刚度增大时，球床壁面附近的流线出现分离和聚合，滞留现象明显，但随着换料的进行，仍会被缓慢卸出堆芯。

从图6可见，球体刚度取E=9.8×106Pa、9.8×107Pa时球流速度矢量的分布相近，靠近底锥区域的球体排列有规律性，间隔为一个球径，这反映出壁面效应的作用较强；而中心主流区域的球流速度明显增大。球体刚度取E=9.8×105Pa时，底锥壁面附近球流仍然受壁面效应影响，但速度分布上没有明显断层。

综合图3～6的结果与分析，球体单元杨氏模量的取值对于球流模拟结果具有明显影响。杨氏模量的大小会直接影响球流的速度分布和流线分布。E=9.8×105Pa的杨氏模量虽然可以大幅提高计算效率，但针对高温堆球床堆芯中的燃料元件进行球流模拟时，已经无法捕捉本应出现的滞留现象，只能得到整体流流型，计算结果并不可信。

2.2 堆积高度对于球床换料过程的影响

对不同高度的堆芯球床的换料过程进行模拟，分别计算了8万、10万、15万、20万个球体的堆积球床。根据2.1节的结论，本节模拟采用9.8×106Pa的杨氏模量。图7为4种球体数量下，球床下部截面球流轴向速度的分布，可以看出，随着模拟球床球体数量与球床总体高度的增加，靠近壁面处的球流速度逐渐增大，滞留现象消失；球床中心的球流速度逐渐减小，球床整体流速分布趋于均匀，球床呈整体下移的整体流流型。

图7 不同球体数量下球床下部球流轴向速度分布Fig.7 Velocity profiles for pebble flow at a high cross section for beds with different number of pebbles

球床中部的轴向速度分布如图8所示，可以看出，随着球体数量的增加，球床整体运动规律与球床下部基本一致。但球体数量为15万和20万的模拟中，除堆芯壁面对临近球体产生拖曳作用外，球流整体轴向流速非常稳定。球体底锥已经无法对30d高度处的球流产生影响。这是由于球体的重力作用随着上方球体数量的增长而愈发显著，底锥对于球流的影响受到削弱。

图8 不同球体数量下球床中部球流轴向速度分布Fig.8 Velocity profiles for pebble flow at a low cross section for beds with different number of pebbles

2.3 卸料过程中球床整体流指数的变化

为进一步分析球床堆芯中的球体滞留现象，引入整体流指数 (mass flow index, MFI)[18]进行研究，并模拟了球床仅卸球不加球的卸料过程，以获得球床在不同球体数量下的流动规律。此时卸球速率仍控制为8个/s，球体杨氏模量取为9.8×106Pa。堆芯由初始0 s时刻的15万个球体，逐步卸出至7 000 s时的9.4万个球体，覆盖了2.2节流型发生转变的中间态。

引入整体流指数MFI为：

(8)

图9为球床不同高度处MFI随换料时长的变化，可以看出，球床底部滞留现象比上部严重：高度为12d和20d的球床底部，由于贴近球床底锥，MFI很低，滞留现象明显；而在高度为30d和40d的较高处，球流趋近于整体流，滞留现象明显减弱。

图9 球床不同高度MFI的变化Fig.9 Mass flow index for pebble beds at several heights

此外，随着卸料过程的推进，堆芯球床内球体数量逐渐减少，MFI呈现逐步下降的趋势。当卸料进行到5 000 s左右时，球床内球体约为11万个，堆芯高约2.8 m，此时在40d高度处，MFI发生明显转变，堆芯随后整体转变为漏斗流流型。模拟卸料过程的初始状态和最终状态分别与2.2节15万个和10万个球体的模拟结果和分析结论吻合。