AP1000核电机组巨型冷却塔型体优化数值计算

周兰欣， 马少帅， 弓学敏， 孙会亮

（1.华北电力大学 电站设备状态监测与控制教育部重点实验室，保定071003；2.保定华电电力设计研究院，保定071003）

近年来核电建设由沿海向内陆地区扩展，大批内陆核电厂的选址、选型工作正在进行.受内陆地区水源限制，核电厂多采用闭式循环的自然通风冷却塔作为冷源设备.与燃煤电厂相比，核电机组循环水量大，所配备的冷却塔淋水面积多在16 000m2以上，塔高超过165m，此类规模的巨型冷却塔在常规火电机组中至今仍无应用.巨型冷却塔的结构参数直接影响循环水的冷却温降，在设计选型过程中，以往规范标准所推荐的经验值有待进一步论证.为了寻找最优结构参数，需要对塔内流场、温度场进行数值计算.

作为第三代核电技术的AP1000机组具有较高的安全可控性，是未来我国核电建设的主流机型之一.根据内陆中部省份地区的年平均气象参数得知，年平均气温为17℃，当冷却塔循环倍率采用50时，AP1000核电机组的循环水量为180 000t／h，所需淋水面积为20 000m2.针对这一淋水面积的冷却塔，采用CFD三维数值计算软件对双曲线塔体进行了模拟计算，以确定最优结构参数.

在设计施工过程中，双曲线冷却塔外型结构根据其塔筒旋转壳体的母线方程进行尺寸放样定位，如图1所示，将横坐标轴取在喉部，其母线方程为

式中：r为塔筒半径；r0为喉部半径；h为坐标点处到喉部的距离（向上为正，向下为负）；α为决定双曲线曲率的特征值［1］.

图1 塔体结构尺寸Fig.1 Structural parameters of the cooling tower

冷却塔的高度由热力计算得出，原则是塔筒抽力应克服塔内阻力，本文中冷却塔塔高为210m.针对AP1000机组，当淋水面积20 000m2确定时，决定双曲线塔筒外型的参数是母线方程中的曲率特征值α及喉部半径r0.以往相关规范标准建议塔筒高度不高于165m，并给出了塔体相关结构尺寸的经验推荐值［2］.显然，在巨型塔的设计、建造过程中，需要对以往经验值进一步验证.

以往常规火电机组冷却塔在选型设计时母线方程中的曲率特征值α常取经验值0.15～0.17.相关标准中喉部面积与塔底面积比的推荐值为0.30～0.40［2］.笔者采用Fluent三维数值模拟软件对核电站用巨型冷却塔母线方程中的曲率特征值α、喉部半径r0的取值重新进行计算验证.α取值为0.14、0.15、0.16、0.17和0.18共5种情况，喉部半径与塔底半径比r0／r2取值为0.5、0.6、0.7和0.8共4种情况.

1 计算方法

计算流体动力学（CFD）是建立在经典流体力学与数值计算方法基础上的一门新兴独立学科，Fluent作为大型商业化CFD软件，目前被广泛应用于数值计算领域.该软件提供了大量计算流体相关模型和算法，并设置了外接用户端口（UDF），具有稳定性好、适用范围广、精度高和拓展性好等优点［3］.因此，选用Fluent作为冷却塔数值模拟平台.

在湿式自然通风冷却塔三维数值计算中，通常分喷溅区、填料区和雨区3个区域对空气与冷却水之间的传热传质过程进行模拟.由于冷却塔内水和空气的体积比小于10%，在喷溅区和雨区采用离散相模型计算塔内流场［4］，其中水滴作为离散相采用拉格朗日法计算，空气作为连续相采用欧拉法求解.填料区采用自定义函数添加源项的方法进行模拟［5－11］，过程是根据填料层的热力特性及气水热质传递原理计算出冷却水的蒸发量和温度降低值，相关计算过程可参阅文献［1］，再通过自定义源项的方式加入到连续相的主控方程中，实现对填料层气水热质传递的耦合计算［4］.

Fluent求解器采用分离隐式算法，流场计算采用典型的Simple算法.计算中，采用稳态雷诺应力平均N－S方程，选用标准k－ε湍流模型，在输运方程中考虑了浮力项.

1.1 湿空气控制方程

冷却塔内湿空气作为连续相，当机组在稳定工况下运行时，其内外流场可以当作稳态计算，通用控制方程如下［4］

式中：ρ为空气密度；ui为速度矢量；Φ为通用变量，分别表示各向速度分量（u、v、w）、湿空气组分Y、温度T、湍动能k和湍流耗散率ε；ΓΦ为广义扩散系数；SΦ为广义源项.

1.2 水滴控制方程

冷却塔内循环水以水滴形式在喷溅区和雨区自由下落，可采用拉格朗日方法计算其流场.离散相模型可追踪水滴的运动轨迹，并耦合计算水滴与湿空气间的热质传递.水滴温度变化方程［4］：

式中：Mp、Tp、Tadb、Ap分别为控制单元内水滴质量、水滴温度、气相干球温度和水滴面积；cp为水滴的比定压热容；h、hfg和t分别为传热系数、水的汽化潜热和经历时间.

水滴蒸发速率为：

式中：cs为水滴表面蒸汽的物质的量浓度；c∞为湿空气中水蒸气的物质的量浓度；Mw为水的摩尔质量；hD为以浓度差为动力的传质系数.

h、hD由 Nusslet关系式得出［4］：

式中：Re为雷诺数；Prva为湿空气的普朗特数；Sc为施密特数.

当Fluent耦合计算气水间热质传递时，相间能量、质量和动量传递将作为源项添加到主控方程中.

1.3 阻力计算

热态的冷却塔阻力包括气水两相间作用力和结构阻力，其中相间作用力对连续相压降的影响通过离散相模型实现，结构阻力主要包括进风口阻力、填料阻力、配水管网阻力和收水器阻力等4部分.由于这些细密的结构体在Gambit中建模困难，因此它们对连续相的作用作为附加动量汇添加到主控方程中，通过用户自定义函数UDF编程实现.

动量汇方程一般形式如下［4］：

式中：Vp为空气通过几何边界面的垂直速度分量；K为压力损失系数.

K由基于试验的经验关系式计算得到［11］

式中：mw为水滴质量流量；ma为湿空气质量流量；L为填料层高度.

Kjfk、Kgw、Kssq分别为进风口、管网和收水器的压力损失系数，根据有关试验资料可以取定值［12］，Kjfk＝0.5、Kgw＝0.5和Kssq＝3.5.

1.4 计算边界及网格划分

以AP1000机组为例，冷却塔高210m，进风口高度12.5m，淋水面积20 000m2，喉部半径与塔底半径比依次选用0.5、0.6、0.7、0.8等4种情况，双曲线塔体母线方程中的特征值α依次选用0.14、0.15、0.16、0.17和0.18等5种情况.塔筒外型依据其母线方程确定，计算环境区域为高度500m、直径500m的圆柱体.

图2给出了计算边界条件.静风条件下进出口均为压力边界条件；塔壁、地面为无滑移壁面边界条件；离散相边界在进出口、水池面、地面均为逃逸边界条件.冷却水在喷淋层采用面射流模型引入冷却塔内，水滴当量直径为3mm，对水滴下落至水池面的运动轨迹进行追踪.

图2 计算边界条件Fig.2 Computational boundary conditions

网格划分如图3所示，利用尺寸函数工具Sizefunction，采用非结构化四面体网格划分，并对网格无关性进行验证，计算了140万、160万、180万网格总数下的平均出塔水温.160万、180万网格数量下的塔内流场无明显变化，所得平均出塔水温相近，考虑到计算时间，选用网格总数为162万.

图3 冷却塔网格划分Fig.3 Grid division for the cooling tower

2 计算结果及分析

2.1 冷却塔内流场分析

经计算，得到三维冷却塔x－y截面上的速度等值线分布图.喉部半径与塔底半径比r0／r2＝0.6、曲率特征值α＝0.16时的速度等值线分布如图4所示，静风环境下塔内流场呈轴对称分布.空气在塔内吸热后温度升高，密度小于塔外空气密度，这种密度差使塔内产生抽吸力.冷却塔外气流在抽吸力作用下经进风口流入塔内，受进风口阻力、填料阻力和雨区阻力的作用，风速逐渐降低.

图5给出了不同喉部半径与塔底半径比下的填料层风速值.由图5可以看出，随着喉部半径与塔底半径比的增大，塔内风速也逐渐增大，当r0／r2＝0.7时，风速达到最大值1.53m／s，当r0／r2超过0.7后，流场发生恶化，风速降低.过大的喉部半径与塔底半径比使气流通过面积增加，速度下降，动压降低，不利于发挥双曲线塔筒对气流的抽吸作用，故r0／r2＝0.7为一临界值.

图4 当r0／r2＝0.6、α＝0.16时，x－y 截面上的速度等值线分布图Fig.4 Contours of air velocity in the x－y plane at r0／r2＝0.6，α＝0.16

图5 填料层风速随喉部半径与塔底半径比的变化Fig.5 Air velocity in fill region vs.r0／r2

2.2 喉部半径与塔底半径比对进塔风量和蒸发水量的影响

图6 进塔风量随喉部半径与塔底半径比的变化Fig.6 Inlet air flow vs.r0／r2

图6为不同喉部半径与塔底半径比下进塔风量的统计值.由图6可以看出，进塔风量随喉部半径与塔底半径比的增大而增大，这是因为当淋水面积确定时，喉部尺寸的增大使冷却塔风筒的容积增大，空气通过量相应增加.当r0／r2＝0.7时，冷却塔内气水比为0.82，此时填料区的容积质散系数达到最大值，由水向空气的传热传质最充分，蒸发水量也达到最大值.图7给出了冷却塔内蒸发水量随喉部半径与塔底半径比的变化，质量传递的趋势可从图中看出.

图7 蒸发水量随喉部半径与塔底半径比的变化Fig.7 Water evaporation rate vs.r0／r2

2.3 塔内温度场和平均出塔水温分析

图8为不同喉部半径与塔底半径比下水池面的温度分布图.由图8可以看出，中心区域温度最高，外围温度较低.随着喉部半径与塔底半径比的增大，中心区域及外围的温度均有所下降，在r0／r2＝0.7时温度达到最低值，此后继续增大喉部半径与塔底半径比，水池面温度反而有所回升.

图8 水池面温度分布（单位：K）Fig.8 Contours of water temperature on pool surface（unit：K）

表1为相同运行工况下冷却塔采用不同喉部半径与塔底半径比时计算所得的平均出塔水温.由表1可以看出，平均出塔水温随喉部半径与塔底半径比的增大呈现先降低后升高的趋势.空气作为冷源对塔内较高温度的喷淋水进行冷却，由2.2节分析可知塔内温度的回升是由于进塔风量的减小造成的.水滴下落时间统计结果如图9所示，水滴下落时间越长，其与空气间的换热就越充分，平均出塔水温就越低.

表1 不同喉部半径与塔底半径比下的平均出塔水温Tab.1 Average outlet water temperature vs.r0／r2

图9 水滴下落时间随喉部半径与塔底半径比的变化Fig.9 Drop time of droplet vs.r0／r2

2.4 曲率特征值α对平均出塔水温的影响

表2为r0／r2＝0.6、相同运行工况下，塔筒母线方程选用不同的曲率特征值α时计算所得的平均出塔水温，可以看出随α的增大，平均出塔水温降低.

表2 r0／r2＝0.6时，不同α下的平均出塔水温Tab.2 Outlet water temperature vs.αat r0／r2＝0.6℃

2.5 综合分析α及r0／r2对平均出塔水温和进塔风量的影响

图10为不同α和r0／r2组合方案下计算所得的平均出塔水温.由图10可以看出，当r0／r2一定时，随特征值α的增大，4种喉部半径取值下的平均出塔水温均呈现下降趋势，且均在α＝0.18处取得最低值.

图10 综合α及r0／r2对平均出塔水温的影响Fig.10 Comprehensive influence ofαand r0／r2on average outlet water temperature

表3给出了不同α及r0／r2组合下计算所得的进塔风量.由表3可知，进塔风量的变化趋势与出塔水温的变化趋势相一致，在α＝0.18处达到最大值.同时可以看出，在计算范围内的任意特征值α下，r0／r2＝0.7均为最佳喉部半径尺寸，此时进塔风量最大，出塔水温最低.当r0／r2大于0.7时，出塔水温开始回升.

表3 综合α及r0／r2对进塔风量的影响Tab.3 Comprehensive influence ofαand r0／r2on inlet air flow kg／s

3 结 论

（1）进塔风量随喉部半径与塔底半径比的增大而增大，当r0／r2＝0.7时进塔风量达到最大值，此时平均出塔水温最低，与r0／r2＝0.5时相比，平均出塔水温降低了0.53K，当r0／r2大于0.7时，平均出塔水温又开始回升.

（2）当r0／r2确定时，平均出塔水温随曲率特征值α的增大而降低，α＝0.18时的平均出塔水温比α＝0.14时的平均出塔水温低0.18K.

（3）在选型设计时，AP1000核电机组用冷却塔的喉部半径与塔底半径比r0／r2应取为0.6～0.7，塔筒母线方程中的曲率特征值α应取0.17～0.18.

［1］赵振国.冷却塔［M］.2版.北京：中国水利水电出版社，2001：36－37.

［2］建设部.GB／T 50102—2003工业循环水冷却设计规范［S］.北京：中国计划出版社，2003.

［3］韩占忠.流体工程仿真计算实例与应用［M］.北京：北京理工大学出版社，2004.

［4］Fluent Inc.User＇s guide［M］.Lebanon，USA：Fluent Inc，2003.

［5］张志刚，王玮，曾德良，等.冷却塔出塔水温的迭代计算方法［J］.动力工程学报，2010，30（5）：372－377.ZHANG Zhigang，WANG Wei，ZENG Deliang，et al.The iterative calculation of cooling tower outlet water temperature［J］.Journal of Chinese Society of Power Engineering，2010，30（5）：372－377.

［6］赵元宾，孙奉仲，王凯，等.十字墙湿式冷却塔冷却特性的数值研究［J］.中国电机工程学报，2009，29（8）：6－13.ZHAO Yuanbin，SUN Fengzhong，WANG Kai，et al.Numerical analysis of cooling performance of wet cooling tower with cross wall［J］.Proceedings of the CSSEE，2009，29（8）：6－13.

［7］POPPE M，ROGENER H.Berechnung von Ruckkuhlwerken［J］.VDI－Warmeatlas，1991（1）：1－15.

［8］WILLIAMSON N，ARMFIELD S，BEHNIA M.Numerical simulation of flow in a natural draft wet cooling tower－the effect of radial thermo fluid fields［J］.Applied Thermal Engineering，2008，28（2）：178－189.

［9］AL－WAKED R，BEHINA M.CFD simulation of wet cooling tower［J］.Applied Thermal Engineering，2006，26（4）：382－395.

［10］KLIMANEK A.Solution of heat and mass transfer in counter flow wet－cooling tower fills［J］.International Communications in Heat and Mass Transfer，2009，36（6）：547－553.

［11］KLOPPERS J C，KROGER D G.Loss coefficient correlation for wet cooling tower fills［J］.Applied Thermal Engineering，2003，23（17）：2201－2211.

［12］周兰欣，蒋波，陈素敏.自然通风湿式冷却塔热力特性数值模拟［J］.水利学报，2009，40（2）：208－213.ZHOU Lanxin，JIANG Bo，CHEN Sumin.Numerical simulation for thermal performance of natural draft wet cooling tower［J］.Journal of Hydraulic Engineering，2009，40（2）：208－213.