过冷沸腾自然对流两相CFD模拟及应用

李松蔚，Riccardo PURAGLIESI，杨 帆，余红星，沈才芬

(1.中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都 610213；2.Paul Scherrer Institute, Villigen 5232, Switzerland)

非能动安全系统有助于简化核电站设计及潜在提高其经济性与安全性。然而，考虑到基于自然循环的弱驱动力，必须采用细致的设计和分析来确保系统执行其预期功能。为此，2004—2008年，国际原子能机构开展了关于“非能动系统的自然循环现象，建模和可靠性”的合作研究项目(CRP)，OECD/NEA的文件[1-2]中对事故期间非能动系统可能发生的大量热工水力现象进行了归类。

大水池中水的行为是关注重点之一，包括温度分层、自然/强迫对流及循环、蒸汽冷凝、上界面传热传质等。实验重点研究温度分层过程及其对传热能力的影响[3]，包括在瑞士保罗谢尔研究所(PSI)PANDA台架上开展的以ESBWR(经济简化沸水堆)为原型的旨在为系统程序及CFD程序验证提供基准的ISP-42实验[4-5]、德国于利希研究中心(FZ Jülich)NOKO台架上开展的冷凝器带热能力综合实验[6-7]、德国亥姆霍兹德勒斯登罗森多夫研究中心(FZDR)TOPFLOW台架上开展的可获得二次侧三维传热特性的倾斜管冷凝实验[8]、美国普渡大学PUMA台架上开展的研究直接接触冷凝及热分层实验[9]等。此外德国FZDR还开展了侧壁加热圆柱水箱实验[10]，通过热电偶及探针测得桶中水的局部温度及空泡份额。

由于水池体积较大，局部区域的传热并不意味着池中温度均匀，三维对流影响传热过程，导致温度分层，根据需要可能需采用计算流体力学(CFD)方法。目前计算仍以系统程序为主，采用CFD的计算较少，Krepper[11]针对NOKO台架实验进行了单相CFD模拟，针对TOPFLOW的两相实验进行了两相CFD模拟[12]，获得两维流场和温场，但未与实验数据进行比较。而ISP-42实验虽适用于CFD计算验证，但基准题计算提交均为系统程序计算结果，并未提交CFD模拟结果[5]。实验测量数据用于TRACE及GOTHIC程序模拟对比，研究结果发现，水池侧的传热计算与实验有显著差异[13]，这是因为TRACE是一维集总参数系统程序，且并非基于池式沸腾开发，而GOTHIC程序具有三维计算能力，但采用非常粗的计算网格，对局部沸腾特性模拟能力有限。因此，采用两相CFD进行模拟比较显得尤为必要。本文基于STAR-CCM+软件，采用两相泡状流及壁面沸腾模型，构建过冷沸腾自然对流两相模拟方法，结合侧壁加热圆柱水箱实验数据对过冷沸腾自然对流进行研究，完成过冷沸腾自然对流两相CFD模拟，应用于ISP-42实验自然循环瞬态模拟研究。

1 过冷沸腾自然对流两相CFD模拟

1.1 侧壁加热圆柱水箱实验

为构建适用的过冷沸腾自然对流两相模拟方法，首先采用德国FZDR实验室的圆柱形水箱过冷沸腾自然对流实验[10，14]进行两相CFD模拟研究与验证。实验在常压下进行，实验装置为一圆柱形水箱，直径0.25 m，高0.25 m，水位高0.21 m。壁面均匀加热。实验中采用热电偶及电导探针测量局部温度及空泡份额，图1为水箱及温度、空泡份额测量位置示意图，通过测量壁面温度，可计算得到实际由壁面到水的随时间变化的加热功率曲线[14]。实验持续时间为1 500 s，数百s后，出现明显的温度分层，1 100 s后，空泡探针探测到气泡产生，在整个实验过程中，水箱中心未监测到气泡产生。

图1 温度及空泡份额测量位置Fig.1 Locations of temperature and volume fraction measurements

1.2 几何建模与网格划分

模拟采用两维轴对称几何建模。几何建模高0.21 m，半径0.125 m。采用结构化网格进行划分，并对加热壁面(右侧)进行网格加密。为避免网格相关性，进行了多种网格划分(表1)。图2示出4种网格方案计算所得的速度，选取t=1 000 s、轴向高度0.1 m处进行径向方向速度对比。图中x轴0.125 m表示靠近加热壁面这一侧，而0.0 m表示靠近圆柱水箱中心轴这一侧。由于仅在近壁面处速度较明显，因此图中仅展示了x=0.110～0.125 m这一区域。由图可知，A方案和B方案计算结果重合，因此选取网格划分方案B可满足网格敏感性要求。图2中计算采用的网格划分轴向节点数80，径向节点数45。其中第1层网格为2.7×10-4m，y+为3，这是为了能成功捕捉相变过程及获取更精细的速度场。

表1 网格划分方案Table 1 Mesh setting case

图2 不同网格方案近壁面液相速度计算结果Fig.2 Calculated liquid velocity in near wall region for different mesh setting cases

1.3 数值模型与不确定性分析

两相CFD过冷沸腾自然循环是瞬态计算，模拟求解基于STAR-CCM+的欧拉-欧拉法[15]。壁面沸腾求解基于Kurul与Podowski[16]提出的壁面热流分配(RPI)模型，在该模型中，来自壁面的给定热通量被分为单相传热、基于成核点和淬火建立的气泡生长蒸发传热、气泡成长到一定大小离开壁面时由水代替原来气泡占据位置所带来的传热影响。模型计算了每部分的热通量和壁温。

为模拟主流区的两相流，液相采用了可实现的k-ε双层湍流模型，双层全y+壁面处理。曳力采用标准阻力系数的Schiller-Naumann关系式，该模型适用于球形固体颗粒、液滴和小直径(球形)气泡。未来将考虑升力和湍流分散力。

采用Ranz和Marshall模型考虑相间传热传质。该模型认为液体温度低于饱和温度则冷凝，液体温度高于饱和温度则蒸发。蒸汽气泡直径设置为单分散尺寸分布，它们在流动区域中的某个位置处的直径，线性依赖于该位置处的液体温度。这种简化处理可能会影响相间界面面积，并因此影响冷凝速率。

由于初始水为静止状态，而此后过冷沸腾产生气泡会逸出或冷凝，因此几何建模左侧设置为轴向对称边界条件，右侧设置为加热壁面，底部为壁面无滑移。水箱上表面设置为气相渗透，该边界条件是使用无滑移壁面及气相渗透壁面对气体逸出表面进行建模，可与自由出口或压力出口边界条件共同使用。STAR-CCM+自动计算通过气相渗透边界的气体质量流量(逸气率)。采用该边界条件，允许气相自由逸出，但不允许液相逸出，考虑到此时处于过冷沸腾状态，水箱中心液体温度仍低于饱和温度，逸出的蒸汽较少，不足以降低自由液面，计算未考虑自由液面高度降低。

自然对流由流体密度不同所驱动，液相须考虑随温度变化导致的密度变化以及浮力作用，因此水物性(密度、比定压热容、动力黏度、热传导率)采用随温度变化的函数，气相密度采用随当地压力变化的函数，浮力采用Boussinesq浮力模型。

由于在低压情况下由压降引起的饱和温度下降较明显，水箱内不同高度压力不同，因此饱和温度被设置为与压力相关的函数。参考压力设置为0.1 MPa。

计算主要关注温度变化，针对与此相关的过冷沸腾模型中的气泡脱离直径进行了敏感性分析，计算对比了模型默认气泡脱离直径、2倍模型默认气泡脱离直径及1/2倍模型默认气泡脱离直径，温度计算结果无明显改变，后继计算模型参数均采用默认值。

从模型研发上来说，采用CFD计算自然循环需考虑对高雷诺数流动可能超出Boussinesq近似的有效范围；雷诺应力和雷诺通量均基于各向同性开发，而自然循环重力起重要作用，需考虑各向异性[3]。计算未出现高雷诺数流动情形，同时已考虑重力影响，模型适用于当前研究情形。

1.4 计算结果与分析

1) 两相CFD模拟结果

图3示出了两相CFD计算不同时刻的水温云图与液相速度场。两相CFD模拟可成功再现实验中由于气泡产生流动加速所导致的第1次温度阶跃(800～1 200 s)，以及由于水箱整体水温达到饱和出现的第2次阶跃(约1 300 s)。模拟展示出随时间的推移，沸腾在加热面水箱上表面出现，由于水箱上部蒸汽产生的影响，自然循环显著，温度分布相同；而下部由于还是单相，仍保持稳定的热分层状态。随时间的推移，温度分布相同的区域和热分层的区域水平界线将逐渐向下推移。两相CFD模拟展示了沸腾出现后水箱中显著的自然循环，如采用单相CFD模拟，由于不考虑相变，无法正确模拟过冷沸腾出现后的温度特性，水箱中不会出现类似自然循环结果，仅在加热面水箱上表面流动较显著。

2) 不同CFD软件计算对比

将本文STAR-CCM+软件计算结果与文献[14]中采用CFX软件计算结果进行对比。图4示出了CFX软件和STAR-CCM+软件典型位置计算温度与实验测量温度随时间的变化对比。从图4可见，CFX对水箱上表面的计算更准确，但对水箱下部的计算误差较大，出现温度第1次阶跃后，CFX的预测值较实验值偏低。STAR-CCM+对水箱上表面温度计算过高的原因是在STAR-CCM+中沸腾模型与壁面滑移无法同时使用，因此计算的水箱上表面速度非常小，这使得温度交混很弱，导致温度被高估。而在实验中，由于水箱上表面为自由表面，设置为可滑移边界条件是更符合实际的。该问题可通过未来STAR-CCM+中在边界构造特殊的速度函数解决。两种软件均具有自然对流过冷沸腾两相模拟能力，STAR-CCM+除水箱上表面温度计算过高外，其余部分温度模拟均比较准确，因此可采用STAR-CCM+进行下一步研究。

图3 两相CFD计算水温云图与液相速度场Fig.3 Two-phase CFD calculation temperature and liquid velocity fields

图4 不同位置计算及测量温度随时间的变化Fig.4 Calculation and experiment temperatures of different locations vs. time

3) 实验验证

用STAR-CCM+对实验进行模拟验证。图5示出了水箱近壁面测量位置计算的温度与实验测量温度随时间变化的对比。从图5可知，计算和实验结果符合良好，合理再现了实验中的两次温度阶跃。

图6示出了计算空泡份额与实验空泡份额随时间变化的比较结果。空泡份额上升趋势及出现位置、出现时间均与实验符合较好。刚出现两相时，仅在接近上表面的加热壁面附近出现气泡，随时间的推移，蒸汽存在区域由壁面向中心迁移。由于实验中空泡探针探头直径为1 mm，无法探测1.5 mm直径以下的气泡，因此计算空泡份额较实验空泡份额高是合理的。

图5 各测量位置计算水温与实验水温的对比Fig.5 Comparison of calculation and experiment water temperatures at different measurement locations

图6 各测量位置计算空泡份额与实验空泡份额的对比Fig.6 Comparison of calculation and experiment void fractions at different measurement locations

通过对圆柱形水箱两相自然对流实验进行两相CFD数值模拟计算，完成了网格敏感性分析、时间步长选择、模型选择，并利用实验数据对数值计算结果进行了验证和评价，数值计算结果与实验数据在温度、空泡出现时间、出现位置上均符合较好，该套计算方法可用于非能动安全壳冷却系统(PCCS)计算。

2 ISP-42实验模拟应用

在PANDA台架上开展的ISP-42实验原型为ESBWR非能动安全壳冷却系统。实验包括6个阶段(A～F)，代表一系列操作模式或过程，每个阶段实际上是1个独立的实验，具有自己的初始条件和边界条件。模拟所关注的是B阶段，由于蒸汽进入管道而导致的PCCS水箱传热沸腾。装置(图7)为高5 m的水箱，中间有20根通蒸汽的管道，水箱通过自然循环带走管道的热量。两相模拟主要关注3 000～3 500 s时间段。在这个时间段，流体温度上升，壁面已有气泡产生，但主流仍无气泡存在[5]。

图7 实验装置流体温度测量位置Fig.7 Liquid temperature measurement location of experiment facility

计算采用二维几何建模，建模选择中间剖面。图8示出了实验中不同监测点水温随时间的变化与计算结果的对比，计算与实验结果符合较好，可较好跟随流体温度变化。再次证明了该两相模拟方法的适用性。蒸汽的产生加速了水的流动，单相时自然对流不强，而两相的出现加强了自然对流。

图8 测量点计算与实验水温随时间的变化Fig.8 Calculation and experiment water temperatures at measured point vs. time

3 结论

通过模拟侧壁加热圆柱形水箱过冷沸腾自然循环实验，成功构建了两相CFD模拟方法，实现了过冷沸腾自然对流瞬态模拟。

1) 两相CFD模拟可获得与实验较一致的温度分布、气泡产生时间与产生位置。可正确模拟沸腾出现后的温度阶跃，以及水箱内显著的自然循环，单相CFD无法模拟该现象。两相CFD计算结果能合理解释实验现象的发生，证实了在过冷沸腾自然对流上进行数值计算的合理性与准确性。

2) 采用相同的两相CFD模拟对PANDA台架上ISP-42实验进行模拟(基于欧洲ESBWR PCCS原型)，可获得与实验较一致的温度分布，再次验证了该模拟方法的正确性。

该两相CFD模拟方法可用于构建自然对流传热特性模型，为非能动安全壳冷却系统及非能动余热排出系统设计提供新的设计手段。