COPRA双层熔池湍流模型优选

Optimization of turbulent models for a corium pool research apparatus （COPRA） double⁃layer melt pool

摘要

为了预测下封头内双层熔池的流动和传热过程，基于不同湍流模型，同时采用凝固熔化模型对堆芯熔池研究装置（corium pool research apparatus，COPRA）双层熔池实验进行计算流体力学（CFD）数值模拟，通过数值计算获得准稳态下熔池的温度、沿壁面的热流密度与内壁面壳层的分布，将模拟结果与实验值进行比较，评价不同湍流模型的适用性和准确性，并进行湍流模型优选。结果表明，壁面模化大涡模拟（WMLES）湍流模型对下封头内双层熔融池流动与传热模拟的准确性和适用性最好；基于WMLES湍流模型，氧化层温度随着熔池高度增大而增大，氧化层上部存在强烈的湍流，在熔池底部的壳层最厚。

Abstract

In order to predict the flow and heat transfer process of a double-layer molten pool at the lower head， computational fluid dynamics （CFD） numerical simulation of a corium pool research apparatus （COPRA） double-layer molten pool experiment has been carried out for different turbulence models by adopting the solidification melting model. Through numerical calculation， the quasi-steady state temperature of the molten pool， the heat flux along the wall and the distribution of the crust were obtained. The simulation results were compared with the experimental values to evaluate the applicability and accuracy of different turbulence models， and the turbulence model was optimized. The results show that the wall-modelled large-eddy simulations （WMLES） turbulence model has the best accuracy and applicability for the simulation of the flow and heat transfer in the double-layer molten pool at the lower head. The WMLES turbulence model shows that the temperature of the oxide layer increases with the increase of the height of the molten pool. Intense turbulence is observed in the upper region of the oxide layer， and the thickest shell forms at the bottom of the molten pool.

Graphical abstract

关键词

双层熔池 / 湍流模型 / 大涡模拟 / 数值模拟

Key words

double-layer melt pool / turbulence models / large-eddy simulation / numerical simulation

引 言

在核电厂发生严重事故后，堆芯以及各类构件由于无法及时冷却发生熔化，堆芯熔融物跌落到压力容器底部，在反应堆压力容器（reactor pressure vessel，RPV）下封头内形成高温熔融池。在堆芯熔融物的作用下，下封头底部由于受到持续强烈的加热可能会发生烧蚀现象。采用外部反应堆压力容器冷却（ERVC）使熔融物在容器内滞留（IVR）是缓解反应堆中发生严重事故的手段之一。

国内外研究者们开展了多种实验来模拟研究压力容器下封头内熔池的自然对流换热行为，如COPO系列实验^［1-2］、BALI实验^［3］、LIVE实验^［4］以及COPRA实验^［5-7］。Zhang等^［5-7］采用1/4圆形切片，按1∶1的比例还原反应堆压力容器下封头，模拟了下封头内熔融物的高瑞利数自然对流行为。Wei等^［8］基于堆芯熔池研究装置（corium pool research apparatus，COPRA）水实验，对不同湍流模型下熔池内自然对流传热特性的变化规律进行数值模拟，并提出将全浮力效应引入到相关的湍流模型中，以更好地模拟自然对流现象。Zhang等^［9］采用壁面模化大涡模拟（WMLES）湍流模型，对BALI、LIVE和COPRA实验进行数值模拟，获得了与实验相吻合的结果，表明WMLES湍流模型能够较好地描述高瑞利数下的自然对流传热行为。由于反应堆中轻金属与氧化物并不相溶，金属层会浮于氧化层之上形成双层熔池结构。Zhang等^［10］利用大涡模拟方法对HPR1000反应堆的单层和双层结构进行了研究，发现双层结构下的壁面熔化比单层结构更为严重。基于COPRA装置，Zhou等^［11］利用10 mm厚的不锈钢板将熔池隔开，研究分层效应对熔池传热特性的影响，以更好地了解双层熔池内的传热规律，然而由于两层采用了同一种模拟材料，不能真实地反映反应堆内的流动情况，故主要侧重于对不同因素的定性研究。目前采用计算流体力学（computational fluid dynamics，CFD）模拟的研究大多针对单层COPRA实验，关于双层COPRA实验的仿真研究较少，因此本文通过数值模拟方法对双层COPRA实验中熔池的流动传热行为进行预测，并基于实验数据对采用的湍流模型进行评价。

1 数学模型

CFD是运用计算机技术研究复杂流动的一种方法，本文将采用Fluent软件进行数值模拟，分析COPRA双层堆芯池自然对流换热实验。基于先前的模拟^［8］，湍流模型的选择对有内热源熔池的自然对流行为具有显著影响。RANS湍流模型在很多研究中得到了广泛应用，因此本文选用RNG k⁃ε模型、Realizable k⁃ε模型、SST k⁃ω模型3种雷诺时均Navier-Stokes（Reynolds average Navier-Stokes， RANS）湍流模型，同时选用尺度自适应模拟（scale adaptive simulation，SAS）模型与WMLES模型进行数值模拟。

1.1 RANS湍流模型

1.1.1 RNG k⁃ε模型

RNG k⁃ε模型在形式上类似于标准k⁃ε模型，但在ε方程中增加了一个附加项，考虑了旋转效应，提高了高速涡流的计算精度，并且RNG k⁃ε模型对近壁区进行适当处理后可以计算低雷诺数效应。总的来说，RNG k⁃ε模型对于强涡流或带有近壁流动的流场模拟具有优越性。

1.1.2 Realizable k⁃ε模型

Realizable k⁃ε模型采用了新的湍流黏度公式，此外ε方程是根据涡度波动均方差的精确运输方程推导出来的。该模型满足雷诺应力的某些数学约束，与真实湍流一致，在模拟自由流动和边界层流动等数值模拟中取得了较好的效果。

1.1.3 SST k⁃ω模型

SST k⁃ω模型在湍流黏度定义中考虑了湍流剪切应力传输，添加涡黏性限制器来解决涡黏性高估的问题。SST k⁃ω模型比标准k⁃ω、BSL k⁃ω（Baseline k⁃ω）模型在壁面附近具有更广的适用范围和更高的准确度。

1.2 SAS湍流模型

SAS湍流模型是RANS湍流模型的改进，具有在不稳定流动条件下捕获湍流能谱的优点，并可对不同频率的涡流进行辨识。由于它在定常流动区域表现出RANS模型的特性，在脉动显著的区域表现出LES模型的特性，因此通常将其归类为RANS/LES混合模型。

1.3 WMLES湍流模型

Zhang等^［10］基于Boussinesq假设描述了LES方法的控制方程，方程式如下。

\frac{\partial \overline{u_i}}{\partial x_i}=\mathrm{0}（1）

\frac{\partial \overline{u_i}}{\partial t}+\frac{\partial \overline{u_i}\overline{u_j}}{\partial x_j}=-\frac{\mathrm{1}}{\rho }\frac{\partial \overline{p}}{\partial x_i}-\frac{\partial {\tau }_{ij}}{\partial x_j}+\frac{\partial }{\partial x_j}·

\left[\mu \left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)\right]-g_i\beta \overline{\varphi }（2）

\frac{\partial \overline{\varphi }}{\partial t}+\frac{\partial {\overline{u}}_j\overline{\varphi }}{\partial x_j}=-\frac{\partial h_j}{\partial x_j}+\frac{\partial }{\partial x_j}\left(\alpha \frac{\partial \overline{\varphi }}{\partial x_j}\right)（3）

{\tau }_{ij}=\overline{u_i{u}_j}-\overline{u_i}\overline{u_j}（4）

q_j=\overline{u_j\varphi }-{\overline{u}}_j\overline{\varphi }（5）

式中，\rho为流体密度，t为时间，x为坐标分量，\overline{u_i}表示平均速度分量，\mu为动力黏度，{\tau }_{ij}为亚网格尺度应力，\overline{p}为平均压强，\beta为膨胀系数，\alpha为热扩散率，\overline{\varphi }为亚网格尺度平均湍流通量，q_j为亚网格尺度通量。亚网格比例参数由以下公式得到。

{\tau }_{ij}-\frac{\mathrm{1}}{\mathrm{3}}{\tau }_{kk}{\delta }_{ij}=-\mathrm{2}{\mu }_t{\overline{S}}_{ij}（6）

{\overline{S}}_{ij}=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)（7）

式中，{\delta }_{ij}为Kronecker三角张量，{\mu }_t为亚网格尺度湍流黏度，{\overline{S}}_{ij}为分辨率尺度的平均应变张量。

WMLES模型能够避免相关层流边界层低雷诺数问题中的雷诺数缩放限制^［12］。涡流黏度可通过混合长度计算。

{\nu }_t=\mathrm{m}\mathrm{i}\mathrm{n}\left[{\left(\kappa d_{\mathrm{w}}\right)}^{\mathrm{2}},{\left(C_{\mathrm{S}\mathrm{m}\mathrm{a}\mathrm{g}}\Delta \right)}^{\mathrm{2}}\right]\cdot S\cdot \left(\mathrm{1}-{\mathrm{e}}^{{}^{-{\left(\frac{y^{+}}{\mathrm{25}}\right)}^{\mathrm{3}}}}\right)（8）

\Delta =\mathrm{m}\mathrm{i}\mathrm{n}\left[\mathrm{m}\mathrm{a}\mathrm{x}\left(C_{\mathrm{w}}{d}_{\mathrm{w}},C_{\mathrm{w}}{h}_{\mathrm{m}\mathrm{a}\mathrm{x}},h_{\mathrm{w}\mathrm{n}}\right),h_{\mathrm{m}\mathrm{a}\mathrm{x}}\right]（9）

式中，d_{\mathrm{w}}为网格中心与壁面间距，参数\kappa=0.418 7，C_{\mathrm{S}\mathrm{m}\mathrm{a}\mathrm{g}}=0.2，C_{\mathrm{w}}=0.15，S是应变率，y^{+}为垂直壁面的无量纲距离；h_{\mathrm{m}\mathrm{a}\mathrm{x}}和h_{\mathrm{w}\mathrm{n}}分别为壁面网格的最大尺度和沿壁面法线网格间距。

1.4 凝固熔化模型

考虑到内壁面壳层的形成，在Fluent的凝固熔化模型中没有直接模拟固-液两相界面，而是将多孔区域方法应用于固液两相糊状区的模拟。熔融材料的液体分数被用作多孔区域中的孔隙率。凝固模型中的能量方程如下。

\frac{\partial }{\partial t}\left(\rho H\right)+\nabla \cdot \left(\rho \mathit{u}H\right)=\nabla \cdot \left(k\nabla T\right)+s（10）

式中，H=h+\mathrm{\Delta }H，h=h_{\mathrm{r}\mathrm{e}\mathrm{f}}+{\int }_{T_{ref}}^T{c}_p\mathrm{d}T，H为任意时刻的焓，h为显热焓，\mathrm{\Delta }H为相变潜热焓，\mathrm{\Delta }H=\beta L，\beta为液相率，L为物质的相变潜热，h_{\mathrm{r}\mathrm{e}\mathrm{f}}为基准焓，c_p为定压比热，k为导热系数，s为源项。

液体分数可以通过以下方程计算。

（11）

如果液体分数大于0且小于1，则将该区域视为糊状区域，动量方程如下。

\frac{\partial }{\partial t}\left(\rho \mathit{u}\right)+\nabla \cdot \left(\rho \mathit{u}\mathit{u}\right)=-\nabla p+\nabla \cdot \left(\tau \right)+\rho g+\mathit{F}（12）

式中，\tau为应力张量。

在动量计算方面，糊状区被视为多孔介质。力项\mathit{F}可以表示为S_m，由方程（13）计算。

S_m=\frac{{\left(\mathrm{1}-\beta \right)}^{\mathrm{2}}}{\left({\beta }^{\mathrm{3}}+\epsilon \right)}{A}_{\mathrm{m}\mathrm{u}\mathrm{s}\mathrm{h}}\left(\mathit{u}-{\mathit{u}}_{\mathrm{p}}\right)（13）

式中，A_{\mathrm{m}\mathrm{u}\mathrm{s}\mathrm{h}}为糊状区常数，{\mathit{u}}_{\mathrm{p}}为固体速度，即凝固的物质通常被连续地从计算域拉出的速度，故也称作拉出速度。

湍流方程中添加的下沉项表示如下。

S_t=\frac{{\left(\mathrm{1}-\beta \right)}^{\mathrm{2}}}{{\beta }^{\mathrm{3}}+\epsilon }{A}_{\mathrm{m}\mathrm{u}\mathrm{s}\mathrm{h}}\phi（14）

式中，\phi 表示要求解的湍流量。由于较高的糊状区常数可能导致模拟过程中出现振荡解，故本文模拟取A_{\mathrm{m}\mathrm{u}\mathrm{s}\mathrm{h}}=3 000。

2 物理模型与数值模拟

2.1 几何模型和边界条件

根据文献［11］中的实验建立双层熔池模型，其氧化层与金属层用隔板隔开，几何模型与边界条件如图1所示。Zhang等^［13］对不同实验中的传热努塞尔数进行比较，发现二维和三维实验结果具有合理的一致性，二者误差均在25%以内，故本文模拟仅考虑二维情况，该方法在单层模拟中也有应用^［8］。容器外壁面为等温冷却，温度为300 K；顶部边界为辐射传热，外部辐射率为0.45，外部辐射温度为实验环境温度300 K，其他边界采用绝热边界。模拟选择双层COPRA实验中的工况1，即氧化层加热功率为10 kW，但由于实验中存在的各种热损失，实际内部加热功率仅为总加热功率的47%，因此氧化层的能量源项为8 030 W/m³。模拟双层均选用20%（摩尔分数，下同）NaNO₃–80%KNO₃的二元混合物，密度随温度的变化基于Boussinesq假设，初始温度为623 K，与实验中盐的浇注温度相同。623 K时材料的相关物性见表1。

2.2 建模及网格划分

整个计算域在ANSYS ICEM中进行网格划分。为了提高网格质量和网格精度，对圆弧几何结构进行L型剖分，并且对壁面处的网格作加密处理，壁面网格最小尺寸为1 mm，使得y⁺＜1，以有效描述壁面处的流体传热行为，如图2所示。

2.3 无关性验证与数值模拟

为了避免计算资源的浪费，选择运算时间较短的RNG k⁃\mathrm{\epsilon }模型，计算时间为1 000 s。以模型的中轴线温度为参考温度，网络无关性验证结果如图3所示。研究发现，网格数对氧化层温度场的影响很大，但对金属层温度场的影响较小。随着网格加密程度增加，氧化层处的温度差异减小。最终采用1.5万网格数模型进行后续的工作。

如图4所示，分别选取1 s、0.2 s和0.02 s这3种时间步长进行步长独立性验证。当步长为0.2 s和0.02 s时，二者的温差不大于1 K，结果基本一致。若步长取得太小会大大增加计算周期，综合考虑，在本文计算中选取时间步长为0.2 s。

本文使用SIMPLE算法对熔池模型进行模拟计算。动量方程、能量方程为二阶迎风离散格式，对于具有高瑞利数自然对流的流动，压力插值使用PRESTO！，总时间为15 h。

3 模拟结果与分析

3.1 模型验证

图5为不同高度下5种湍流模型的熔池温度与实验值的比较。由于实验采用热电偶测量，故温度是同一高度上所有热电偶的平均值。可以看出，5种模型的最大误差都出现在底部区域，且WMLES模型的误差最小，这是由于目前的几何模型没有考虑壳层和下封头之间存在的气体间隙。

在同一高度下模拟值和实验值存在误差，这是由于在实验中氧化池内采用电加热棒加热，加热区域不均匀，而数值模拟假设整个氧化池为均匀衰变热；其次COPRA双层实验采用三维几何形状，而数值模拟采用二维几何形状。此外，其他影响因素还包括数值误差、热电偶的不准确性等。然而尽管所有湍流模型与实验都存在着一定的误差，但总体趋势一致。

定义无量纲参数Θ和Φ以作进一步分析。Θ表示某一处模型值与实验值的拟合度。

\Theta =\left|\frac{X_{\mathrm{c}\mathrm{a}\mathrm{l}}-X_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}{X_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}\right|（15）

式中，X_{\mathrm{c}\mathrm{a}\mathrm{l}}和X_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}分别为相同状态下的模型值和实验值。

Φ表示模型值与实验值的整体拟合度。

\Phi =\sqrt[]{\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{\Theta }_i^{\mathrm{2}}}（16）

式中，n为监测值数量。

Φ越趋近于0，代表模型值与实验值拟合度越高。表2列出了不同湍流模型下温度的{\Theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}值、{\Theta }_{\mathrm{m}\mathrm{i}\mathrm{n}}值、{\Theta }_{\mathrm{a}\mathrm{v}\mathrm{e}}值和\Phi值，可以看出，WMLES模型的拟合度最高。

为了更好地描述传热特性，定义无量纲参数q/q_{\mathrm{a}\mathrm{v}\mathrm{e}}来对比数值模拟结果与实验结果，其中q为每个极角下的局部热流密度，q_{\mathrm{a}\mathrm{v}\mathrm{e}}为沿弯曲壁面的平均热流密度。q/q_{\mathrm{a}\mathrm{v}\mathrm{e}}沿极角的分布如图6所示。在30°以下时，壁面热流密度很小且几乎维持不变，这是因为在熔池底部内表面形成壳层，降低了热流密度，从而有利于IVR的实现。随着角度的增加，热流密度迅速上升，最大热流密度出现在不锈钢板下面（图6绿色虚线）。但是在整个金属层中，由于顶部的辐射换热，局部热流密度有所减小，降低了壁面的热流密度。

由图6可以看出，5种模型的总体趋势与实验值一致，主要差异在氧化层的上部，其中RNG k⁃ε、Realizable k⁃ε和WMLES模拟值的最大热流密度与实验值更接近。由于采用的熔盐热导率远低于实际金属材料的热导率，实验和模拟中均没有发现金属层的聚焦效应以及下封头的烧蚀。

表3列出了不同湍流模型下壁面热流密度的{\Theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}值、{\Theta }_{\mathrm{m}\mathrm{i}\mathrm{n}}值和\Phi值。可以看出，模拟获得的最大热流密度显著大于实验值，WMLES模型的\Phi值均最小，拟合度最高，这与Zhang等^［9］对单层熔池模拟的结果相同。

对极角60°以下壳层厚度的模拟结果和实验结果进行对比，见图7。随着角度的增加，壳层厚度逐渐变薄。模拟值与实验值存在一定的误差，这是由于实验中加热棒的长度相对较短，导致弯曲壁面附近加热不足，从而产生了较厚的壳层。在模拟中，最厚的壳层形成于熔池底部，约0.3 m，而在实验中约0.23 m，但模拟与实验结果的总体趋势是一致的。另外，在30°以下区域壳层较厚，所以在低角度处热流密度呈现出平坦的分布，这与图5的热流密度分布一致。但随着角度的增大，壳层厚度逐渐减小，在50°时几乎消失，表明在氧化层上部区域存在着较强的湍流以及很高的壁面热流密度。

表4列出了不同湍流模型下壳层厚度的{\Theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}、{\Theta }_{\mathrm{m}\mathrm{i}\mathrm{n}}、{\Theta }_{\mathrm{a}\mathrm{v}\mathrm{e}}和\Phi值。可以看出，RNG k⁃ε和Realizable k⁃ε模型有着相似的{\Theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}值、{\Theta }_{\mathrm{m}\mathrm{i}\mathrm{n}}值，而SST k⁃ω和SAS模型相似，尽管WMLES模型有着最小的{\Theta }_{\mathrm{a}\mathrm{v}\mathrm{e}}值，但RNG k⁃ε模型的整体拟合度最高。

对比5种模型的计算结果，综合考虑温度、热流密度和壳层厚度的拟合度情况，选择WMLES模型进行下封头内双层熔池的数值模拟。

3.2 温度场

图8显示了计算时长为1 h、3 h、6 h、9 h、12 h和15 h的双层熔池温度场。温度场可以分为两个部分，即底层的热分层和中上部的均匀温度分布。熔融混合物在容器壁的冷却及重力作用下，沿着内壁面流动，并在熔池底部形成涡流，将底部的热量置换到上部，所以氧化层底层的温度低于中上部的温度。随着池高的增加，熔池温度逐渐增加，但是由于顶部的辐射冷却，金属层的温度低于氧化层。整个熔池中，在同一高度下壁面处的温度更低。由于金属层不产热，随着容器壁冷却与顶部的辐射冷却，整个金属层会发生熔盐凝固。从图8中可以看出，12 h和15 h的温度场相似且温度差在1 ℃以内，可认为15 h时刻双层熔池已经达到准稳态，故下文仅研究15 h时刻准稳态下的速度场与液相场。

如图9所示，由于底部区域湍流较弱，在容器壁冷却的情况下形成壳层，其热阻较大，导致了更大的温度梯度。随着池高的增加，温度上升速度降低，温度场趋于平缓，熔池温度在不锈钢板下达到峰值（图9绿色虚线）。然而，由于金属层不产生热量，温度在金属层中迅速下降，在准稳态条件下，温度下降到与环境温度一致。

3.3 速度场

图10显示了15 h时刻双层熔池速度场，图11为双层熔池速度矢量图。对于氧化层来说，速度场可以分为3个区域，即上部湍流混合区，沿容器壁的向下流动区和下部静态区。熔盐沿着容器壁向下流动，由于驱动力不足，少量的熔盐在壁面处形成涡流，而剩余部分仍流向熔池底部。氧化层上部区域存在不规则的流动行为，湍流强烈，故温度较为均匀。对于靠近壁的区域，流动速度先增加后减小，在涡流周围存在着热分层。而下部区域的熔盐速度非常小，且有横向流动趋势，从传热角度看，该区域主要以热传导为主，这是该区域存在热分层现象的主要原因。

3.4 液相场

图12显示了15 h时刻氧化层的液相体积分数分布。在重力和自然对流的共同作用下，凝固的熔盐向熔池底部迁移，壳层厚度随着熔池壁面径向角度增大而减小。氧化层中下部区域均处于糊状区，仅上部还依然保持纯液态状态。由图10的速度分布可知，在约50°时，壳层的尖端附近产生涡流，使该区域的温度分布更加均匀。另外，壳层还会影响容器内壁面的热流密度，充当容器和熔池之间的热阻，有效阻止下封头壁面发生熔化，因此，氧化层底部的热流密度小于顶部的热流密度。

4 结论

（1）在Fluent平台上基于COPRA双层实验，对双层熔池自然对流传热过程进行数值模拟。在温度和热流密度上，WMLES模型的\Phi值均最小，拟合度最高；而在壳层厚度上，尽管WMLES模型有着最小的{\Theta }_{\mathrm{a}\mathrm{v}\mathrm{e}}值，但RNG k⁃ε模型的整体拟合度最高。综合考虑实验与模拟结果，选择WMLES模型进行下封头内双层熔池的数值模拟计算，误差最小。

（2）稳态下在氧化层下部发生热分层现象，中上部温度均匀分布，熔池温度随着熔池高度增大而增大，但金属层的温度比氧化层低，整个金属层会发生熔盐凝固现象。

（3）氧化层上部存在强烈的湍流，壳层尖端附近存在小部分的涡流，而下部区域的熔盐流动速度非常低，且具有横向流动趋势。

（4）壳层厚度随着极角角度的增加而逐渐减小，在熔池底部壳层最厚。

参考文献

[1] KYMÄLÄINEN O， TUOMISTO H， HONGISTO O， et al. Heat flux distribution from a volumetrically heated pool with high Rayleigh number［J］. Nuclear Engineering and Design， 1994， 149（1-3）： 401-408.

[2] KYMÄLÄINEN O， TUOMISTO H， THEOFANOUS T G. In‑vessel retention of corium at the Loviisa plant ［J］. Nuclear Engineering and Design， 1997， 169（1-3）： 109-130.

[3] BONNET J M. Thermal hydraulic phenomena in corium pools： the BALI experiment［C］//Proceedings of the Workshop on Severe Accident Research： Workshop on Severe Accident Research Held in Japan （SARJ-98）. Tokyo， 1999： 79-86.

[4] BUCK M， BÜRGER M， MIASSOEDOV A， et al. The LIVE program – results of test L1 and joint analyses on transient molten pool thermal hydraulics［J］. Progress in Nuclear Energy， 2010， 52（1）： 46-60.

[5] ZHANG L T， ZHANG Y P， ZHOU Y K， et al. COPRA experiments on natural convection heat transfer in a volumetrically heated slice pool with high Rayleigh numbers［J］. Annals of Nuclear Energy， 2016， 87： 81-88.

[6] ZHANG Y P， ZHANG L T， ZHOU Y K， et al. The COPRA experiments on the in⁃vessel melt pool behavior in the RPV lower head［J］. Annals of Nuclear Energy， 2016， 89： 19-27.

[7] ZHANG Y P， ZHANG L T， ZHOU Y K， et al. Natural convection heat transfer test for in⁃vessel retention at prototypic Rayleigh numbers – results of COPRA experiments［J］. Progress in Nuclear Energy， 2016， 86： 80-86.

[8] WEI H Y， CHEN Y T. Assessment of different turbulence models on the large scale internal heated water pool natural convection simulation［J］. Annals of Nuclear Energy， 2019， 131： 23-38.

[9] ZHANG L T， LUO S M， ZHANG Y P， et al. Large eddy simulation on turbulent heat transfer in reactor vessel lower head corium pools［J］. Annals of Nuclear Energy， 2018， 111： 293-302.

[10] ZHANG L T， ZHOU Y K， LUO S M， et al. Large eddy simulation for the thermal behavior of one⁃layer and two⁃layer corium pool configurations in HPR1000 reactor［J］. Applied Thermal Engineering， 2018， 145： 38-47.

[11] ZHOU Y K， WU S H， ZHANG Y P， et al. Experimental research on heat transfer behavior of large scale two⁃layer salt melt pool based on COPRA facility［J］. Annals of Nuclear Energy， 2020， 138： 107166.

[12] SHUR M L， SPALART P R， STRELETS M K， et al. A hybrid RANS-LES approach with delayed⁃DES and wall⁃modelled LES capabilities［J］. International Journal of Heat and Fluid Flow， 2008， 29（6）： 1638-1649.

[13] ZHANG L T， ZHOU Y K， ZHANG Y P， et al. Natural convection heat transfer in corium pools： a review work of experimental studies［J］. Progress in Nuclear Energy， 2015， 79： 167-181.