一种获得偏离泡核沸腾型临界热流密度的方法
【专利摘要】一种获得偏离泡核沸腾型临界热流密度的方法，步骤如下：1.设置一个初始热流密度qm，为了获得用于临界热流密度机理模型所需的局部热工水力参数，基于均匀流模型建立管内两相热工水力瞬态求解模型，基于有限差分法来求解两相均匀流偏导，以及采用半隐式差分格式求解离散后的控制方程；2.当步骤1计算所得局部热工水力参数达到稳态时，将步骤1获得的局部热工水力参数作为已知参数赋给偏离泡核沸腾型临界热流密度机理模型，求解临界热流密度qCHF；3.判断qm与qCHF是否一致，若没有收敛，则更新qm的值，重复步骤1，直至二者相同；本发明方法具有较好的拓展性，对数据依赖性低，能够有效地预测临界热流密度，对反应堆安全分析具有重要意义。
【专利说明】
一种获得偏离泡核沸腾型临界热流密度的方法
技术领域
[0001] 本发明属于反应堆安全分析技术领域，具体涉及一种获得偏离泡核沸腾型临界热 流密度的方法。 技术背景
[0002] 为保证堆芯流通通道的正常运行和完整性，需要对管内流动换热特性进行研究， 临界热流密度是反应堆设计的主要热工安全准则之一。临界热流密度是核燃料元件表面发 生传热恶化时的热通量，是冷却剂流动沸腾机理发生转变的结果。如果反应堆燃料元件表 面发生沸腾临界，将会导致燃料元件表面温度过高从而造成加热壁面烧毁，放射性物质泄 露，造成严重的运行事故，直接影响到反应堆的安全运行。因此准确预测临界热流密度对反 应堆的设计和安全分析、对反应堆系统综合性能的提高具有重要意义。
[0003] 目前，临界热流密度的预测方法较多采用适用于不同流体和不同结构流道的各种 经验关系式以及临界热流密度查询表。经验关系式是最早使用的预测和研究临界热流密度 的一种方法。至今，针对不同的热工水力条件、流型、管道结构尺寸等参数，国内外研究者们 已经提出数百种临界热流密度预测经验关系式。但大多数经验关系式是基于一个比较小的 数据库开发得到，因此大多数经验关系式适用范围有限，不能扩展至适用范围以外的工况。 而临界热流密度查询表可在很宽的参数范围直接查表并插值的方法得到相应参数下的临 界热流密度，但其只适用于单一流体，且外推到查询表提供的参数范围之外(例如高流速和 高过冷度)时预测精度下降，较经验关系式需要花费更多的计算时间。此外，数据表格的形 式使其不便于更新。
[0004] 基于此，有必要提供一种可避免使用大量的实验数据拟合公式，降低对数据的依 赖性，并能够适用于不同的流型和流体工质的更有效率的预测临界热流密度的方法。
[0005] 发明说明
[0006] 为了解决上述问题，本发明提供一种获得偏离泡核沸腾型临界热流密度的方法， 可避免使用大量的实验数据拟合公式，降低对数据的依赖性，并能够适用于不同的流型和 流体工质。
[0007] 为了达到上述目的，本发明的采用如下技术方案：
[0008] -种获得偏离泡核沸腾型临界热流密度的方法，包括如下步骤：
[0009] 步骤一:给定流道结构尺寸以及边界条件(进口温度、进口流量和出口压力），设置 一个初始热流密度qm，采用均匀流模型模拟流道内两相瞬态流场，以此获得流通通道中的 局部热工水力参数，即压力、流量及焓值，从而为偏离泡核沸腾型临界热流密度计算提供参 数；
[0010] 步骤二：当步骤一中计算所得局部热工水力参数达到稳态时，将步骤一得到的出 口处局部热工水力参数作为已知参数赋给偏离泡核沸腾型临界热流密度机理模型，采用该 机理模型计算得到临界热流密度qcHF;;
[0011]步骤三：比较步骤一中的初始热流密度qm与步骤二计算所得临界热流密度qcHF，， 若qm大于qcHF，则减小qm，若qm小于qcHF，则增加 qm，采用变化后的热流密度重复步骤一的过 程;若二者相同，则得到所求工况下的临界热流密度qcHF;
[0012] 步骤一所述的采用均匀流模型模拟流道内两相瞬态流场，采用有限差分法求解两 相均匀流偏导，采用半隐式差分格式求解离散方程。
[0013] 步骤二所述的偏离泡核沸腾型临界热流密度机理模型为微液层蒸干模型，其原理 为:在壁面热流密度逐渐缓慢增加过程中，流体达到两相，壁面附近会出现一个拉长的汽 块，它由多个小汽泡结合而成;长汽块下裹覆着一层很薄的液体层(微液层），该长汽块会沿 加热壁面移动，此时长汽块会阻碍主流液体对壁面的有效冷却；当长汽块下的微液层蒸干 时，该点处的传热会剧烈恶化，从而导致壁温急剧升高，进而发生沸腾临界;求解机理模型 的关键公式如下：
[0014]
[0015] 式中：
[0016] Pf--饱和液相密度/kg · m-3;
[0017] hfg--汽化潜热/J .kg-1;
[0018] δ一一长汽块下微液层厚度/m;
[0019] Ub--长汽块移动速度/m · s-1;
[0020] Lb--长汽块长度/m。
[0021] 从式（1)中可知，微液层厚度(δ)、汽块移动速度(Ub)和汽块长度(Lb)是求解的关键 参数；
[0022 ]其中汽块移动速度Ub通过求解汽块的轴向受力平衡来确定，轴向方向上汽块的受 力包括浮力Fb和拖曳力Fd，公式如下：
[0023] Fb-Fd = O 公式（2)
[0024]微液层厚度δ通过求解汽块的径向受力平衡来确定，径向方向上的汽块受力平衡 包括侧面提升力Fr、蒸发力Fe3以及壁面润滑力Fw:
[0025] -Fr+Fe+Fw=0 公式（3)
[0026]汽块长度Lb采用Helmholtz临界波长，其公式如下所示：
[0027]
[0028] 式中：
[0029] P1--液相密度/kg · m-3;
[0030] pg一一汽相密度/kg· m-3;
[0031 ] σ--表面张力/N · m工。
[0032] 本发明和现有技术相比，具有如下优点：
[0033] 1.采用均匀流模型模拟了两相瞬态流场，可较精确地计算模拟出临界沸腾发生点 处的热工水力参数。
[0034] 2.从微观上对汽泡受力平衡和运动进行研究，得到了较为精细的机理模型，同时 也便于后续分析运动条件对模型的影响因素。
【附图说明】
[0035] 图1为本发明方法的流程图。
[0036] 图2为一维空间网格。
[0037]图3为汽块受力的示意图。
[0038]图4为机理模型计算流程图。
【具体实施方式】
[0039]下面结合附图和【具体实施方式】对本发明结构进行详细说明。
[0040]如图1所示，本发明一种获得偏离泡核沸腾型临界热流密度的方法，具体方法如 下：
[0041 ]步骤一:给定流道结构尺寸以及边界条件(进口温度、进口流量和出口压力），设置 一个初始热流密度qm，采用均匀流模型模拟流道内两相瞬态流场，以此获得流通通道中的 局部热工水力参数，即压力、流量及焓值，从而为偏离泡核沸腾型临界热流密度计算提供参 数；
[0042] 步骤二：当步骤一中计算所得局部热工水力参数达到稳态时，将步骤一得到的出 口处局部热工水力参数作为已知参数赋给偏离泡核沸腾型临界热流密度机理模型模块，计 算得到临界热流密度qcHF;
[0043]步骤三：比较步骤一中的初始热流密度qm与步骤二计算所得临界热流密度qcHF，若 qm大于qcHF，则减小qm，若qm小于qcHF，则增加 qm，采用变化后的热流密度重复步骤一的过程； 若二者相同，则得到所求工况下的临界热流密度qcHF。
[0044] 下面进行详细说明：
[0045] 步骤一中：两相瞬态流场采用均匀流模型模拟求解，下面给出具体的模型：
[0046] 由于所研究的临界热流密度特性只需要提供外部的热工水力参数，对流场模拟的 要求不高，因此采用改进的均匀流模型模拟流场。均匀流模型的基本假设为：
[0047]①两相具有相等的线速度，即Ug = Ul = UH;
[0048]②两相之间处于热力平衡；
[0049]③使用合理确定的单相摩擦因数。
[0050]运用上述假定于两相流动基本方程，并与单相流体守恒方程类比，定义合适的混 合物物性和摩擦因数计算式，可以得到以下守恒方程。
[0051 ] 1)质量守恒方程：
[0052]
公式（5)
[0053] 式中：
[0054] p一一两相混合物的密度/kg · Hf3;
[0055] u一一两相混合物的流动速度/m · ?Γ1;
[0056] t--时间/s;
[0057] z--长度/m。
[0058] 考虑两相流体的可压缩性，为了求解方程组，增加如下物性方程：
[0059] p = p(p,h)公式（6)
[0060]由式(6)求导可得：
[0061]
[0062]
[0063]
[0064]
[0065]
[0066]
[0067]
[0068]
[0069]
[0070]
[0071]
[0072] 式中：
[0073] Prh 加热周长/m;
[0074] A--流通面积/m2;
[0075] qv--体积释热率/W · m 3;
[0076] ①若无内热源，即流体内部没有释热。上式右侧的qv为零。
[0077] ②忽略动能变化引起的能量变化，则上式左侦
[0078]因此，能量方程可简化为
[0079]
[0080]
[0081]
[0082]
[0083]
[0084]步骤一中：上述守恒方程，采用有限差分法求解偏导，采用半隐式差分格式求解离 散方程。
[0085]在空间上，控制体划分采用交错网格。交错网格中，热力学变量(压力p，焓值h，密 度P，空泡份额α等)存储在控制体中心（图2中虚线位置），而流体的速度存放在控制体边界 处(图2中实线位置）。
[0086]在时间上采用半隐式差分格式进行离散。为了加快计算速度，除了连续方程和能 量方程中对流项的速度，以及动量方程中的压力梯度项采用隐式格式以外，其他项均采用 显示格式。
[0087] 各守恒方程的离散形式如下：
[0088] 席醫Φ栢卞'fS.
[0089]
[0090]
[0091]
[0092]
[0093]
[0094]
[0095] 时间步长的选取可以违反(u+a) Δ t/ Δ Z<1的波速限制，但是要满足u Δ t/ Δ Z<1 的流速限制。
[0096] 能量守恒方程：
[0097]
[0098]
[0099]
[0100] 其中Prh 加热周长/m。
[0101] 交错网格以及有限差分法的优点是压力梯度项与速度项自然的耦合在一起，并且 如果满足Courant法则u Δ t/ Δ Z<1则形成一个对角优势矩阵，可保证求解算法的稳定性。 所有网格的方程最后可以组成一个压力的三对角矩阵，它可以利用直接消去法求解，也可 以利用Gauss-Seidel迭代法求解。总体计算思路是:①用动量方程求解出速度;②利用上一 时层求解压力方程;③用解出的压力回代求解速度和焓值。
[0102] 步骤二中：微液层蒸干机理的具体模型如下：
[0103] 微液层蒸干模型假设加热壁面附近产生的小汽泡结合形成大汽块，在汽块下存在 非常薄的液相层，称为微液层。汽块移动过程中，当汽块下的液体全部蒸发烧干时，该点处 的加热壁面被单相蒸汽覆盖从而导致传热恶化，进而导致沸腾临界发生。求解机理模型的 关键公式如下：
[0104] qcHF = 5pfhfgUb/Lb 公式（19)
[0105] 式中：
[0106] pf--饱和液相密度/kg · m-3;
[0107] hfg一一汽化潜热/J .kg-1;
[0108] δ一一汽块下微液层厚度/m;
[0109] Ub--汽块移动速度/m · s-1;
[0110] Lb--汽块长度/m。
[0111] 从（14)式中可知，微液层厚度(δ)、汽块移动速度(Ub)和汽块长度(Lb)是求解的关 键参数。
[0112]如图3所示，上述三个关键参数可以通过分析汽块受力获得，其中汽块移动速度Ub 通过求解汽块的轴向受力平衡来确定。轴向方向上汽块的受力包括浮力Fb和拖曳力Fd，公式 如下：
[0113]
[0114]
[0115]
[0116]
[0117] 式中：
[0118] Db--汽块当量直径/m;
[0119] P1--液相的密度/kg · m-3;
[0120] pg一一汽相的密度/kg· m-3;
[0121] g一一重力加速度/miT2;
[0122] t--时间/s;
[0123] Cd--拖拽系数；
[0124] Ubi一一汽块中心线所处的径向位置处的主流速度/ms'
[0125] 联立式(20)-(22)，汽块速度Ub可表达如下：
[0126]
[0127]
[0128]
[0129]
[0130]
[0131] 式中：
[0132] σ--表面张力/N · m S
[0133] D--管道直径/m;
[0134] Tw--壁面剪切力/kg · m-1 · s-2。
[0135 ] 采用Karman速度分布方程确定当地流速：
[0136]
[0138] 近壁面区的汽块一般处在缓冲区的范围，即对应于Karman速度分布方程的第二个 关系式。所以为简化计算，可选用第二个关系式计算当地流速分布，该关系式可以改写为：
[0137]
[0139]
[0140]
[0141]
[0142]
[0143]
[0144] 微液层厚度δ通过求解汽块的径向受力平衡来确定。径向方向上的汽块受力包括 平衡侧面提升力Fr、蒸发力Fe3以及壁面润滑力F w:
[0145] -Fr+Fe+Fw=0 公式（30)
[0146] 蒸发力是由微液层中液相蒸发产生的小汽泡进入汽块，对汽块冲击而产生的，这 个力将推动汽块向远离加热壁面。侧面提升力是由汽块和主流区存在速度差异，使汽块在 运动过程中发生旋转而产生，提升力的大小还与汽块在径向方向上的速度梯度有关，它会 推动汽块靠近加热壁面。
[0147] 沸腾临界将要发生时，由微液层蒸发而产生的蒸发力为：
[0148]
公Λ (31)[0149] 垂直管道内两相湍流流动中的侧面提升力的表达式：
[0150；
[0151；
[0152；
[0153；
[0154] 其中式(33)中，α为临界沸腾发生点处的空泡份额，式(29)中的速度由式(22)计算 获得。
[0155] 因此，侧面提升力可表达如下：
[0156]
[0157] 此外还需要考虑壁面润滑力Fw，壁面润滑力将非常靠近壁面处的汽泡推向中心 区，而对离壁面距离较远的汽泡作用力则非常小，其表达式如下：
[0158]
[0159]
[0160]
[0161] 上式中，&和(：2为经验常数，yw表示离加热壁面的垂直距离，式(36)在y w〈5Db范围内 有效。
[0162] 联立式（30)，（31)，（32)和(35)，微液层厚度δ可表达如下：
[0163]
[0164] 求出δ以后，可以通过式(28)计算得到Ubl，然后分别用式(23)和(29)计算得到U b和 U，再将新的Ub和Lb代入式(38)计算得到新的δ，通过迭代计算直到以上几个参数值达到收 敛。最后，将迭代计算得到的S，U b和Lb代入式（19)来计算临界热流密度值。机理模型的流程 图如图4所示。
【主权项】
1. 一种获得偏离泡核沸腾型临界热流密度的方法，其特征在于:包括如下步骤： 步骤一:给定流道结构尺寸以及边界条件即进口温度、进口流量和出口压力，设置一个 初始热流密度qm，采用均匀流模型模拟流道内两相瞬态流场，以此获得流通通道中的局部 热工水力参数，即压力、流量及焓值，从而为偏离泡核沸腾型临界热流密度计算提供参数； 步骤二：当步骤一中计算所得局部热工水力参数达到稳态时，将步骤一得到的出口处 局部热工水力参数作为已知参数赋给偏离泡核沸腾型临界热流密度机理模型，采用该机理 模型计算得到临界热流密度qcHF; 步骤三:比较步骤一中的初始热流密度qm与步骤二计算所得临界热流密度qcHF，，若qm大 于qCHF，则减小qm，若qm小于qCHF，则增加 qm，采用变化后的热流密度重复步骤一的过程;若二 者相同，则得到所求工况下的临界热流密度qcHF。2. 根据权利要求1所述的一种获得偏离泡核沸腾型临界热流密度的方法，其特征在于： 步骤一所述的采用均匀流模型模拟流道内两相瞬态流场，采用有限差分法求解两相均匀流 偏导，采用半隐式差分格式求解离散方程。3. 根据权利要求1所述的一种获得偏离泡核沸腾型临界热流密度的方法，其特征在于： 步骤二所述的偏离泡核沸腾型临界热流密度机理模型为微液层蒸干模型，其原理为:在壁 面热流密度逐渐缓慢增加过程中，流体达到两相，壁面附近会出现一个拉长的汽块，它由多 个小汽泡结合而成；长汽块下裹覆着一层很薄的液体层即微液层，该长汽块会沿加热壁面 移动，此时长汽块会阻碍主流液体对壁面的有效冷却；当长汽块下的微液层蒸干时，该点处 的传热会剧烈恶化，从而导致壁温急剧升高，进而发生沸腾临界;求解机理模型的关键公式 如下： qcHF = 8pfhfgUb/Lb 公式（1) 式中： Pf--饱和液相密度/kg ? nf3; hfg一一汽化潜热/J ? kg' 8一一长汽块下微液层厚度/m; Ub--长汽块移动速度/m ? s S Lb--长汽块长度/m。 从式(1)中可知，微液层厚度S、汽块移动速度Ub和汽块长度Lb是求解的关键参数； 其中汽块移动速度Ub通过求解汽块的轴向受力平衡来确定，轴向方向上汽块的受力包 括浮力Fb和拖曳力Fd，公式如下： Fb-Fd = 0 公式(2) 微液层厚度S通过求解汽块的径向受力平衡来确定，径向方向上的汽块受力平衡包括 侧面提升力Fr、蒸发力以及壁面润滑力Fw: -Fr+Fe+Fw=0 公式(3) 汽块长度Lb采用Helmholtz临界波长，其公式如下所示：Pi一一液相密度/kg ? m-3; Pg一一汽相密度/kg ? nf3; 8--表面张力/N ? nf1。
【文档编号】G06F17/50GK106055850SQ201610566357
【公开日】2016年10月26日
【申请日】2016年7月18日
【发明人】苏光辉, 黄思洋, 张大林, 田文喜, 秋穗正
【申请人】西安交通大学