-
多孔介质模型主要通过在动量方程中附加源项来模拟固体区域对于流动的影响,源项包括粘性损失和惯性损失两部分。在本文的模拟中,计算多孔介质区的压降方程为
$$\Delta p = - \sum\limits_{j = 1}^3 {{C_{{2_{ij}}}}\left( {\frac{1}{2}\rho {v_j}\left| {{v_j}} \right|} \right)} ,$$ (1) 其中:
$\Delta p$ 是多孔介质区的在$j$ 方向的压降;${C_2}$ 是多孔介质区惯性阻力系数;$\rho $ 是液态铅铋的密度;${v_j}$ 是液态铅铋在$j$ 方向的流速。在FLUENT中设置多孔介质的惯性阻力因子,就可以模拟多孔介质区域的压降损失。根据各区域的几何结构与流动状态来估算压降损失,根据多孔介质模拟的压降与之相等,从而得到惯性阻力因子
${C_2}$ 。${C_2}$ 的计算通过FLUENT用户自定义函数(UDF)实现,不同结构采用不同的压降模型。 -
根据各区域的结构特征和相应的流动状态,得到液态铅铋流经各区域的阻力系数计算公式,即为各结构的压降模型,如表1所列。
表 1 堆内压降模型
项目 纵向阻力系数 横向阻力系数 参考文献 配重区、操作头区、换热器 $\begin{aligned}&{f_{\rm{a}}} = 0.048R{e^{ - 0.2}}\\&{C_{\simfont\text{纵}}} = 4\frac{{{f_a}}}{{{D_{\rm{h}}}}}\frac{1}{{{\gamma ^2}}} \end{aligned}$ $\begin{aligned}& Eu{\rm{ = 1}}{\rm{.256\,39}} \times \beta \times R{e^{ - 0.183\,81}}\\& {C_{\simfont\text{横}}} = \frac{{Eu \times n}}{{a \times {\gamma ^2}}}\end{aligned}$ [4-5] 主泵、流量分配器、围筒出口中部多孔板 $\xi = \left\{ \begin{array}{l} 7.68{\left[ { {\beta ^{ - 1} }{ {\left(\dfrac{t}{d}\right)}^{ - 0.035} } - 1} \right]^2}{\beta ^{0.94} }{\rm{ } }\quad{{Re} } < 1.1 \times {10^5} \\ 4.50{\left[ { {\beta ^{ - 1} }{ {\left(\dfrac{t}{d}\right)}^{ - 0.052} } - 1} \right]^2}{\beta ^{0.58} }{\rm{ } }\quad{{Re} } \geqslant 1.1 \times {10^5} \\ \end{array} \right.$ 忽略 [6] 组件入口段与出口段 局部阻力损失的经验或理论公式:
进口段:2×90°弯头+2×沿程+1×多孔板+1×突扩;
出口段:2×117.5°折管+2×沿程+1×多孔板。忽略 [7] 燃料组件组件
含绕丝棒束$\begin{aligned}f = \frac{{64}}{{Re}}{F^{0.5}} + \frac{{0.081\,6}}{{R{e^{0.133}}}}{F^{0.935\,5}}\frac{{{N_r}\pi \left( {D + {D_w}} \right)}}{{{S_t}}}\\F = {\left( {\frac{P}{D}} \right)^{0.5}} + {\left[ {7.6*\frac{{D + {D_w}}}{H}{{\left( {\frac{P}{D}} \right)}^2}} \right]^{2.16}}\end{aligned}$ $\begin{aligned}&Eu{\rm{ = 1}}{\rm{.256\,39}} \times \beta \times R{e^{ - 0.183\,81}}\\&{C_{\simfont\text{横}}} = \frac{{Eu \times n}}{{a \times {\gamma ^2}}}\end{aligned}$ [5, 8-9] 哑组件中段 尼古拉兹实验圆管湍流半经验公式 忽略 [7] -
FLUENT多孔介质模型的传热模型分为热平衡模型和非热平衡模型,二者的差别在于是否将多孔介质的流体与固体的温度视作相等。为了准确地模拟冷却剂与燃料棒束之间的换热,燃料组件使用了非热平衡模型。CiADS铅基堆采用含绕丝的燃料组件,因此选取Pacio等[10-11]对含绕丝组件铅铋流动换热实验的拟合换热系数h作为换热模型:
$$Pe = RePr = \rho {u_b}{d_h}{C_p}{\lambda ^{ - 1}},$$ (2) $$Nu = 2.483 + 0.0312P{e^{0.8}},$$ (3) $$h = \frac{{Nu \times {\lambda}}}{L},$$ (4) 式中:λ为导热系数;L为特征长度;ub为冷却剂主流速度即沿燃料棒束方向的速度;dh为水力直径;Cp为定压比热容。
-
CiADS铅基堆一回路冷却剂为液态铅铋合金,FLUENT计算采用变物性参数,所用的物性关系式如表2所列。
表 2 液态铅铋主要物性关系式
项目 公式 单位 参考文献 密度 ${\rho _{\rm{LBE}}} = 11\,096 - 1.323\,6T$ kg/m3 [12] 导热率 ${\lambda _{\rm{LBE}}} = 3.61 + 0.015\,17T - 0.000\,001\,741{T^2}$ W/(k·m) [12] 定压比热容 ${C_{{P_{\rm{LBE}}}}} = 159 - 0.027\,2T + 0.000\,007\,12{T^2}$ J/(kg·K) [12] 动力粘度 ${\mu _{\rm{LBE}}} = 0.000\,494{{\rm e}^{754.1/T}}$ kg/(m·s) [12] 体膨胀系数 ${\alpha _{\rm{LBE}}} = \dfrac{1}{{8\,383.2 - T}}$ 1/K [13] -
CiADS铅基堆燃料组件活性区在轴向的功率分布近似服从余弦分布[2],采用UDF分段函数模拟热源的空间分布。停堆后衰变热随时间的变化同样通过UDF函数进行描述。
-
本文湍流模型选取标准κ-ε模型,稳态计算的离散格式为Coupled,瞬态计算的离散格式为PISO。由于自然循环为体积力驱动的物理现象,因此压力离散格式选择Body force weighted。此外,为了更准确地求解边界层的流动与传热,近壁面处理采用Enhanced wall treatment。
-
三维模型使用网格数量为84,156,386,655,1 022万的5套多面体网格进行无关性分析。冷却剂循环质量流量与回路的压降密切相关,因此,质量流量达到网格无关性即表明回路的流动特性达到网格无关性,故选取通过某流通面的质量流量作为网格无关性检验的关键参数。该质量流量随网格数量的变化曲线如图5所示,曲线从第三点开始渐渐趋于平稳,第三点与第五点的相对偏差为0.683%,差距较小,因此选择386万的网格。
-
CiADS次临界反应堆额定工况下,堆芯功率为7.74 MW,主换热器和主泵为运行状态,一回路冷却剂循环类型为强迫循环主导的混合循环。额定工况的热工水力设计值及与三维模拟结果如表3所列,相对偏差在0.5%以下,这表明CFD较好地模拟了额定工况。
表 3 三维额定工况计算结果与设计值对比
项目名称 设计值 模拟值 相对误差/% 冷却剂总质量流量 541 kg/s 541.88 kg/s 0.16 哑组件旁通流量
百分比3.88% 3.87% 0.26 燃料组件进口温度 553.15 K 552.43 K 0.07 燃料组件出口温度 653.15 K 652.05 K 0.25 堆芯内冷却剂温升 100 K 99.62 K 0.38 图6图给出了包含主泵截面的速度分布云图,最大流速出现在主换热器与泵的连接管处,为0.61 m/s;泵出口至冷池中部的流速在0.12 ~ 0.18 m/s之间,冷池下部的流速降至0.03 ~ 0.08 m/s。图7为热池的速度矢量图,可以看到换热器两侧的进口存在回流,换热器中形成了两个明显的涡旋;在围筒上端隔板及围筒出口的阻力的影响下,组件出口上方形成了较大的涡旋。图8为包含换热器截面的温度分布云图,冷却剂最高温度出现在堆芯组件上部,为656.3 K;冷池中的冷却剂温度分布均匀,且与热池的温度分界明显,表面围筒与冷热池隔板的热分隔效应良好;燃料包壳的最高温度为662.6 K,对应热源的轴向分布,燃料区呈现出明显的温度梯度。
各结构的静压压降如图9图所示,从图中可见环腔回路的主要压降分布在堆芯组件进出口及棒束区,换热器、流量分配器、围筒等结构的压降损失相比较小。
为了研究额定工况下自然循环对冷却剂循环流量的贡献,设置无重力(无自然循环)为对照组进行计算。无重力额定工况的冷却剂流量为530.96 kg/s,与额定工况相差2.02%。这是由于额定工况下冷却剂流速较大使得回路压降损失也大,从而自然循环驱动压头的相对贡献较低。
-
为了研究CiADS铅基堆在低功率水平运行且主泵故障的情况下,一回路是否能够建立起自然循环保证反应堆安全运行,我们进行了20%额定功率下的三维稳态计算。20%功率运行的主要热工水力参数如表4所列。由于主泵停转,冷却剂循环流量大幅下降,自然循环流量最终稳定在92.53 kg/s,为额定流量的17.10%;流量大幅减少后,流动状态的变化影响了额定工况下完成的堆芯流量分配,旁通流量占比减小了约0.47%,燃料组件之间出口冷却剂的温差相应增大(如表4所列)。
表 4 三维20%功率运行模拟结果
项目名称 模拟值 冷却剂总质量流量/(kg/s) 92.53 哑组件旁通流量百分比/% 3.41 堆芯进口温度/K 507.84 堆芯出口温度/K 605.67 堆芯内冷却剂温升/K 97.83 由表5可见,低功率运行下包壳最高温度与冷却剂最高温度均低于安全限值。即在低功率水平运行时,主泵故障后反应堆仍能依靠一回路建立的自然循环安全运行,表明CiADS铅基堆具有低功率自然循环运行的能力和一定的事故容错能力。
表 5 三维20%功率运行燃料组件温度与质量流速的模拟结果
燃料
组件出口
温度/K冷却剂最高
温度/K固体最高
温度/K质量流速/
(kg/s)FA1 610.22 616.52 619.07 3.38 FA2 612.44 617.38 620.27 3.86 FA3 611.83 617.42 620.30 3.81 FA4 610.45 616.77 619.32 3.35 FA5 608.95 617.11 619.99 3.83 FA6 597.52 614.28 616.30 2.55 FA7 595.77 614.04 615.96 2.42 FA8 596.75 613.80 615.70 2.39 -
由于堆内结构复杂,三维模型网格数量较大,进行瞬态计算的计算量难以接受。因此,本文尝试建立CiADS铅基堆二维等效模型进行瞬态计算。根据体积等效、换热面积等效等原则,对堆内各结构进行等效转换,建立二维旋转对称模型如图10所示。
燃料组件与哑组件的二维等效模型与三维结构的对比如图11所示,堆内各结构的等效原则见表6。
表 6 反应堆二维建模等效原则
项目 等效原则 主容器、换热器 总体积与高度相等、换热面积相等 配重区、操作头区、主泵 高度、体积相等 组件入口段、出口段 高度、流通面积相等 组件中段 高度、流通面积、体积相等 围筒 总体积相等 流量分配器 高度与厚度相等 堆芯支撑板、冷热池隔板 半径与厚度相等 -
对比二维模型网格数分别为50万和80万的两套结构化网格的计算结果,主要热工水力参数如表7所列,可见两套网格的计算结果相差较小,考虑到计算的准确性与经济性,选择50万的网格进行后续计算。
表 7 两套网格计算结果对比
项目名称 50万网格计算结果 80万网格计算结果 冷却剂总质量流量/ kg/s 533.99 534.28 哑组件旁通流量百分比/% 3.96 4.13 燃料组件进口温度/K 553.35 553.09 燃料组件出口温度/K 653.28 652.93 堆芯内冷却剂温升/K 99.93 99.84 -
进行二维等效转换时,很难同时满足所有等效原则,因此有必要对比二维与三维模拟结果,评价二维等效方法带来的偏差。表8~9对比了二者的稳态计算结果,泵动量源项的相对偏差最大,为14.49%,最大绝对偏差为15.21 K,表明二维等效方法具有可行性,基本可以反映三维模拟的结果,因此后续的瞬态分析均采用二维模型进行计算。
表 8 二维与三维额定工况计算结果参数对比
匹配参数 二维 三维 相对偏差 泵动量源项/(N/m3) 301 100 263 000 14.49% 换热器换热系数/(W/m2/K) 930 1 030 9.71% 表 9 二维与三维无重力额定工况与20%功率运行结果对比
项目名称 二维无重力额定工况 三维无重力额定工况 偏差 二维20%功率运行 三维20%功率运行 偏差 冷却剂总质量流量/(kg/s) 534.25 530.96 0.62% 97.66 92.53 5.54% 堆芯进口温度/K 579.91 572.14 7.77 498.84 507.84 9 堆芯出口温度/K 683.93 675.78 8.15 611.88 605.67 6.21 堆芯内冷却剂温升/K 104.02 103.64 0.38 113.04 97.83 15.21 包壳最高温度/K 697.80 686.91 10.89 618.18 620.30 2.12
Study on Natural Circulation and Residual Heat Removal Capability of the Lead-based Fast Reactor
-
摘要: 自然循环特性是铅基反应堆一回路的关键运行特性,对反应堆的非能动应急余热排出具有重要的影响,自然循环特性与余热排出能力是反应堆热工水力研究的重要内容。采用多孔介质方法,建立了CiADS铅基堆1/4三维计算模型,使用FLUENT程序对额定工况与低功率工况进行稳态计算。为了研究全厂断电事故下的余热排出过程,从热工水力的等效原则出发,尝试建立二维等效模型以提高瞬态计算效率。结果表明,CiADS铅基堆具备低功率自然循环运行能力和一定的事故容错能力;二维等效模型与三维模型计算结果吻合较好,可用于瞬态下的简化分析;CiADS铅基堆的非能动余热排出系统能够较好地应对全厂断电事故,反应堆具有良好的固有安全性。Abstract: The natural circulation characteristics is the key operating characteristic of the primary circuit of the lead-based reactor, which has an important influence on the passive residual heat removal system of the reactor. The natural circulation characteristics and residual heat removal capacity in the reactor are an important part of the thermal hydraulic research of reactor. Using the porous media method, the 1/4 3D calculation model of the CiADS lead-based reactor was established, and the steady state calculation of rated conditions and low power operating condition were performed by FLUENT code. In order to study the transient process of residual heat removal of CiADS lead-based reactor under the station blackout accident, a 2D equivalent model which according to the principles of thermal hydraulic equivalence was established for transient calculation. The results show that CiADS lead-based reactor has low power natural circulation operation capability and a certain degree of accident tolerance; The results of 2D equivalent model and 3D model agree well with each other and the 2D model can be used in transient analysis; The 2D transient calculation results show that the passive residual heat removal system can cope well with the station blackout accident, and the CiADS lead-based reactor has good passive safety performance.
-
Key words:
- CiADS /
- lead-based fast reactor /
- natural circulation /
- residual heat removal /
- CFD
-
表 1 堆内压降模型
项目 纵向阻力系数 横向阻力系数 参考文献 配重区、操作头区、换热器 $\begin{aligned}&{f_{\rm{a}}} = 0.048R{e^{ - 0.2}}\\&{C_{\simfont\text{纵}}} = 4\frac{{{f_a}}}{{{D_{\rm{h}}}}}\frac{1}{{{\gamma ^2}}} \end{aligned}$ $\begin{aligned}& Eu{\rm{ = 1}}{\rm{.256\,39}} \times \beta \times R{e^{ - 0.183\,81}}\\& {C_{\simfont\text{横}}} = \frac{{Eu \times n}}{{a \times {\gamma ^2}}}\end{aligned}$ [4-5] 主泵、流量分配器、围筒出口中部多孔板 $\xi = \left\{ \begin{array}{l} 7.68{\left[ { {\beta ^{ - 1} }{ {\left(\dfrac{t}{d}\right)}^{ - 0.035} } - 1} \right]^2}{\beta ^{0.94} }{\rm{ } }\quad{{Re} } < 1.1 \times {10^5} \\ 4.50{\left[ { {\beta ^{ - 1} }{ {\left(\dfrac{t}{d}\right)}^{ - 0.052} } - 1} \right]^2}{\beta ^{0.58} }{\rm{ } }\quad{{Re} } \geqslant 1.1 \times {10^5} \\ \end{array} \right.$ 忽略 [6] 组件入口段与出口段 局部阻力损失的经验或理论公式:
进口段:2×90°弯头+2×沿程+1×多孔板+1×突扩;
出口段:2×117.5°折管+2×沿程+1×多孔板。忽略 [7] 燃料组件组件
含绕丝棒束$\begin{aligned}f = \frac{{64}}{{Re}}{F^{0.5}} + \frac{{0.081\,6}}{{R{e^{0.133}}}}{F^{0.935\,5}}\frac{{{N_r}\pi \left( {D + {D_w}} \right)}}{{{S_t}}}\\F = {\left( {\frac{P}{D}} \right)^{0.5}} + {\left[ {7.6*\frac{{D + {D_w}}}{H}{{\left( {\frac{P}{D}} \right)}^2}} \right]^{2.16}}\end{aligned}$ $\begin{aligned}&Eu{\rm{ = 1}}{\rm{.256\,39}} \times \beta \times R{e^{ - 0.183\,81}}\\&{C_{\simfont\text{横}}} = \frac{{Eu \times n}}{{a \times {\gamma ^2}}}\end{aligned}$ [5, 8-9] 哑组件中段 尼古拉兹实验圆管湍流半经验公式 忽略 [7] 表 2 液态铅铋主要物性关系式
项目 公式 单位 参考文献 密度 ${\rho _{\rm{LBE}}} = 11\,096 - 1.323\,6T$ kg/m3 [12] 导热率 ${\lambda _{\rm{LBE}}} = 3.61 + 0.015\,17T - 0.000\,001\,741{T^2}$ W/(k·m) [12] 定压比热容 ${C_{{P_{\rm{LBE}}}}} = 159 - 0.027\,2T + 0.000\,007\,12{T^2}$ J/(kg·K) [12] 动力粘度 ${\mu _{\rm{LBE}}} = 0.000\,494{{\rm e}^{754.1/T}}$ kg/(m·s) [12] 体膨胀系数 ${\alpha _{\rm{LBE}}} = \dfrac{1}{{8\,383.2 - T}}$ 1/K [13] 表 3 三维额定工况计算结果与设计值对比
项目名称 设计值 模拟值 相对误差/% 冷却剂总质量流量 541 kg/s 541.88 kg/s 0.16 哑组件旁通流量
百分比3.88% 3.87% 0.26 燃料组件进口温度 553.15 K 552.43 K 0.07 燃料组件出口温度 653.15 K 652.05 K 0.25 堆芯内冷却剂温升 100 K 99.62 K 0.38 表 4 三维20%功率运行模拟结果
项目名称 模拟值 冷却剂总质量流量/(kg/s) 92.53 哑组件旁通流量百分比/% 3.41 堆芯进口温度/K 507.84 堆芯出口温度/K 605.67 堆芯内冷却剂温升/K 97.83 表 5 三维20%功率运行燃料组件温度与质量流速的模拟结果
燃料
组件出口
温度/K冷却剂最高
温度/K固体最高
温度/K质量流速/
(kg/s)FA1 610.22 616.52 619.07 3.38 FA2 612.44 617.38 620.27 3.86 FA3 611.83 617.42 620.30 3.81 FA4 610.45 616.77 619.32 3.35 FA5 608.95 617.11 619.99 3.83 FA6 597.52 614.28 616.30 2.55 FA7 595.77 614.04 615.96 2.42 FA8 596.75 613.80 615.70 2.39 表 6 反应堆二维建模等效原则
项目 等效原则 主容器、换热器 总体积与高度相等、换热面积相等 配重区、操作头区、主泵 高度、体积相等 组件入口段、出口段 高度、流通面积相等 组件中段 高度、流通面积、体积相等 围筒 总体积相等 流量分配器 高度与厚度相等 堆芯支撑板、冷热池隔板 半径与厚度相等 表 7 两套网格计算结果对比
项目名称 50万网格计算结果 80万网格计算结果 冷却剂总质量流量/ kg/s 533.99 534.28 哑组件旁通流量百分比/% 3.96 4.13 燃料组件进口温度/K 553.35 553.09 燃料组件出口温度/K 653.28 652.93 堆芯内冷却剂温升/K 99.93 99.84 表 8 二维与三维额定工况计算结果参数对比
匹配参数 二维 三维 相对偏差 泵动量源项/(N/m3) 301 100 263 000 14.49% 换热器换热系数/(W/m2/K) 930 1 030 9.71% 表 9 二维与三维无重力额定工况与20%功率运行结果对比
项目名称 二维无重力额定工况 三维无重力额定工况 偏差 二维20%功率运行 三维20%功率运行 偏差 冷却剂总质量流量/(kg/s) 534.25 530.96 0.62% 97.66 92.53 5.54% 堆芯进口温度/K 579.91 572.14 7.77 498.84 507.84 9 堆芯出口温度/K 683.93 675.78 8.15 611.88 605.67 6.21 堆芯内冷却剂温升/K 104.02 103.64 0.38 113.04 97.83 15.21 包壳最高温度/K 697.80 686.91 10.89 618.18 620.30 2.12 -
[1] WIDER H U, CARLSSON J, LOEWEN E. Progress in Nuclear Energy, 2005, 47(1-4): 44. doi: 10.1016/j.pnucene.2005.05.003 [2] 彭天骥, 顾龙, 王大伟, 等. 原子能科学技术, 2017, 51(12): 2235. doi: 10.7538/yzk.2017.51.12.2235 PENG Tianji, GU Long, WANG Dawei, et al. Atomic Energy Science and Technology, 2017, 51(12): 2235. (in Chinese) doi: 10.7538/yzk.2017.51.12.2235 [3] FLUENT INC. ANSYS 19 Fluent User's Guide.[CP/OL]. Fluent User’s guide, 2019. [4] 黄兴华, 王启杰, 陆震. 化工学报, 2000, 51(03): 297. doi: 10.3321/j.issn:0438-1157.2000.03.003 HUANG Xinghua, WANG Qijie, LU Zhen. Journal of Chemical Industry and Engineering, 2000, 51(03): 297. (in Chinese) doi: 10.3321/j.issn:0438-1157.2000.03.003 [5] 李飞, 孙奉仲, 史月涛, 等. 山东大学学报(工学版), 2014, 44(04): 70. doi: 10.6040/j.issn.1672-3961.0.2014.033 LI Fei, SUN Fengzhong, SHI Yuetao, et al. Journal of Shandong University (Engineering Science), 2014, 44(04): 70. (in Chinese) doi: 10.6040/j.issn.1672-3961.0.2014.033 [6] 赵锴. 静电除尘器内气流均匀性及压力特性试验研究[D]. 杭州: 浙江大学, 2017. ZHAO Kai. Research on Flow Uniformity and Pressure Characteristics of Electrostatic Precipitator[D]. Hangzhou: Zhejiang University, 2017. (in Chinese) [7] 孔珑. 工程流体力学[M]. 北京: 中国电力出版社, 2014. KONG Long. Engineering Fluid Mechanics[M]. Beijing: China Electric Power Press, 2014. (in Chinese) [8] REHME K. Nuclear Technology, 1973, 17(1): 15. doi: 10.13182/NT73-A31250 [9] YAN J, KOCHUNAS B, HURSIN M, et al. Coupled Computational Fluid Dynamics and MOC Neutronic Simulations of Westinghouse PWR Fuel Assemblies with Grid Spacers[C]. The 14th International Topical Meeting on Nuclear Reactor Thermal hydraulics, NURETH-14, 2011. [10] PACIO J, DAUBNER M, FELLMOSER F, et al. Nuclear Engineering and Design, 2016, 301: 111. doi: 10.1016/j.nucengdes.2016.03.003 [11] PACIO J, DAUBNER M, FELLMOSER F, et al. Nuclear Engineering & Design, 2018, 330: 225. doi: 10.1016/j.nucengdes.2018.01.034 [12] KOLOSZAR L, BUCKINGHAM S, PLANQUART P, et al. Nuclear Engineering and Design, 2017, 312: 256. doi: 10.1016/j.nucengdes.2016.05.008 [13] CRESPO L S. Handbook on lead–bismuth Eutectic Alloy and Lead Properties, Materials Compatibility, Thermal-hydraulics and Technologies. Chapter 6: Compatibility of Structural Materials with LBE and Pb: Standardisation of Data Corrosion Mechanism and Rate. OECD/NEA, 2007.