直接空冷凝汽器空气侧热流性能的跨尺度模拟

程通锐，杜小泽，杨立军

(电站设备状态监测与控制教育部重点实验室(华北电力大学), 北京市 102206)



直接空冷凝汽器空气侧热流性能的跨尺度模拟

程通锐，杜小泽，杨立军

(电站设备状态监测与控制教育部重点实验室(华北电力大学), 北京市 102206)

针对冬季直接空冷凝汽器的防冻工业需求，将翅片管和空冷单元这2个尺度的模型耦合起来，既考察翅片管束的整体流动换热性能，又考察直接冷却翅片管内凝结水的冲击换热特性，为确定冬季低温运行条件下空冷凝汽器管内凝结换热甚至结冰计算提供第3类热边界条件。建立了空冷单元模型和翅片管模型，然后，将三维进口风条件下的翅片管级数值模拟结果简化成输入-输出代理模型；从空冷单元系统级数值模拟中获得系统级边界-迎风面的速度分布；进而将翅片管数值模拟代理模型应用到系统边界上，得到了系统边界上翅片管冲击换热特性以及翅片管平均对流换热特性的空间分布。数值计算结果表明，冲击换热系数比平均对流换热系数大一个数量级；在迎风面上冲击换热系数和平均对流换热系数均呈现出上低下高、左右不对称的特点；冲击换热系数随着转速的下降而下降。

直接空冷凝汽器；跨尺度模拟；传热；冻结

Cross-Scale Simulation for Air-side Thermo-Fluid

Performance of Direct Air-Cooling Condenser

0 引 言

直接空冷凝汽器以环境空气作为传热介质冷却汽轮机排汽，具有显著的节水效益，在我国北方富煤缺水地区火力发电中得到广泛应用。但是，空冷火电机组在北方冬季低温环境中运行存在冻结风险。揭示低温环境空气在直接空冷凝汽器单元中的热流行为，是确定空冷机组防冻策略的前提和基础。

空冷单元风机来风以绝对速度对管束进口扁平管圆弧头部内的凝结水冲击冷却，且此处空气接近环境温度，换热强度大，冻结风险高，因而冲击冷却是本文的研究重点。带连续翅片的单排扁平管束是目前大型火电机组直接空冷凝汽器的基本传热元件。有大量文献对扁平管翅片空间空气侧的对流换热进行了研究[1-14]，采用的湍流模型包括标准k-ε模型、低雷诺数模型[10]、kω-SST模型[14]，以及直接大涡模拟方法等[11]。这些模拟结果在不同程度上揭示出空冷翅片管束在不同迎面风速下的传热性能，以及空气流场的流动细节，但从未有文献描述扁平管圆弧部分的冲击换热性能。另外一些研究，则着重于翅片管束及风机组成的空冷凝汽器Λ型单元内流动换热性能[15]，以及阵列运行的轴流风机的集群空气动力学特性。单一的空冷单元级数值模拟只能体现翅片管束的整体性能，难以体现直接冲击冷却管内凝结水的局部冲击换热特性。因而本文将翅片管尺度和空冷单元尺度的数值计算模型结合，同时考察翅片管束的整体换热性能和局部冲击换热性能。

本文以实际结构的直接空冷凝汽器单元为对象，通过建立小尺度翅片管模型，获得翅片管整体性能和关心的小尺度模型细节特征，并采用降维算法将其重构成输入-输出代理模型；建立系统级空冷单元大尺度模型，将小尺度翅片管模型整体性能嵌入大尺度模型进行系统级数值模拟，从而获得系统级边界上的参数分布；然后将关心的小尺度翅片模型的细节特征代理模型应用到系统边界上，从而获得系统边界上的翅片细节特征的分布。

1 物理数学模型

1.1 直接空冷凝汽器单元

本文以600 MW典型空冷凝汽器单元为研究对象，图1显示了1个空冷单元的三维图。左右2组管束成“Λ”型结构布置，夹角为60o，每面管束的整体尺寸为10.8 m×9.81 m，在空冷凝汽器单元的模拟中简化为多孔介质区域。轴流风机位于凝汽器翅片管束正下方。风机直径8.91 m，轴功率80 kW。

假设流动是三维不可压缩稳态流动。为了便于处理浮升力的作用，热物性只考虑密度随温度的变化，守恒量方程中只考虑动量方程中与体积力有关的项，并采用大气温度作为参考温度。湍流模型采用了标准k-ε二方程模型。采用相对坐标系和实体风机叶片相结合的方法建立风机模型。关于空冷单元翅片管束的计算，文献中采用的模型有2种，体散热器(heat exchanger) 模型和散热器平面(radiator)与多孔介质相结合的模型，本文采用了后者。在radiator

图1 直接空冷凝汽器单元物理模型

平面上设置翅片管束的压降与法向速度的关系为

(1)

式中：Δp为翅片管沿程压降；ρ为空气密度；va,n为散热器平面法向风速；kL为无量纲流动阻力系数，表达式为

(2)

式中fn为多项式系数，f1=75.21，f2=-32.5，f3=3.60。散热器平面的热流密度可以表达为

q=h(Ts-Ta,d)

(3)

式中：q为热流密度；h为基于Ta,d定义的等效传热系数；Ts为凝结温度；Ta,d为散热器平面下游相邻网格的温度。等效传热系数的表达式为

(4)

式中hn为多项式系数，h1=537，h2= 2 017，h3=-98。

连续性方程为

(5)

式中：u为空气流速，x为空间坐标。

动量方程为

(6)

式中：μeff为有效粘性系数；p为流体压力；g为重力加速度；S为源项，可用于Radiator条件的表达。

能量方程为

(7)

式中：T为流体温度；e为工质比内能；Sh为源项，可用于Radiator条件的表达。

标准k-ε模型为

(8)

式中：k为湍动能；ε为耗散率；μt为湍流粘性；σk为关于k的湍流普郎特数；Gk表示由速度梯度引起的湍动能的产生项；Gb表示由浮力引起的湍动能的产生项。

控制方程为

(9)

式中：σε为关于ε的湍流普朗特数；C1ε、C2ε、C3ε为模型常数。

如图1所示，空冷单元下方四周围设置为压力入口边界条件。空冷单元四周设置为对称边界条件。空冷单元顶部出口设置为压力出口边界条件。Radiator模型中管束传热系数和流动阻力采用试验关联式[17]，温度采用汽轮机背压对应的饱和温度。对于多孔介质区域，根据通流的方向设置非通流方向的其他两个方向的流动阻力系数为极大值。风机叶片所在的区域设置风机转速等参数。

在计算区域采用有限容积法对控制方程进行离散，对流项采用中心差分，采用了Simple压力速度耦合算法进行迭代求解。为保证计算结果的网格独立性，要求不同网格下计算得到的空气流量变化至稳定，最后采用的模型的总网格数为228万个。在单元级数值模拟中，得到的是空冷单元强迫对流与空气密度差形成的自然对流相结合的混合对流，得到了空气流场以及风机的流量、换热量。

1.2 连续翅片扁平管

风机来风进口条件复杂，因而通过最大转速下单元级的数值模拟得到了翅片管束进口风分速度的最大变化范围(0～10m/s)，而后在此范围内，进行多种分速度搭配下翅片管传热性能的数值模拟。

翅片空气通道的结构如图2所示，基管设为定壁温边界条件，上下面设为周期性边界条件，翅片设置为可与流体换热的耦合壁面条件。

图2 翅片管空气侧计算区域及边界条件

为保证数值模拟进出口处于充分发展的状态，设置了进口段和出口段。进口段四周侧面设置为周期性边界条件。进口段正面设置为速度入口边界条件。空气由于经过了较长的翅片空间的流动诱导，出口流速方向平行于流道，因而出口段四周侧面设置为对称边界条件。出口正面设置为压力出口条件。控制方程及其求解方法类似于单元数值模拟，不同的是湍流模型采用了kω-SST模型，这种模型的特点是可以实现转捩过程的数值模拟，而空冷翅片管的设计风速正处于转捩过渡的范围，且这种模型在文献[14]中得到了实验验证。kω-SST湍流模型方程为

(10)

Gω-Yω+Dω

(11)

最后采用的具有网格无关性的模型的总网格数为105万个。表达流道整体换热性能的传热系数，即平均对流换热系数为

(12)

式中：Q表示换热量；Ar表示基管换热面积；Af表示翅片面积；ηf表示肋化效率；Δtm为对数平均温差。可通过设置翅片温度为基管温度的数值模拟得到理想最大换热量，实际换热量与该理想换热量之比即为肋化效率。表达管束进口段极限换热性能的冲击换热系数，即冲击换热系数为

(13)

式中：qarc,max表示圆弧冲击点的热流密度，同时也是圆弧区域的局部最大热流密度；Tw表示基管温度；T0表示环境温度。

2 跨尺度模拟策略

本文使用Fluent软件对2个尺度的模型及相应耦合策略进行空冷单元空气侧换热的多尺度计算。

(1)获得小尺度模型的整体性能，翅片管的平均对流特性和流动阻力特性已经通过试验获得。

(2)将翅片管的整体性能作为等效子模型嵌入大尺度模型，建立大尺度系统级模型，从而可以进行系统级数值模拟，具体将翅片管平均对流换热特性和流动阻力特性嵌入直接空冷单元模型，获得迎风面的三维风速的分布及其变化范围。

(3)在系统边界条件的变化范围(0～10 m/s)内进行小尺度模型的数值模拟，获得所关心的小尺度模型细节特征，并采用一定算法对关心的细节建立代理模型。即在迎风面的三维风速的变化范围内对翅片管进行数值模拟，获得圆弧头部的冲击换热特性以及作为对比的平均换热特性，然后将冲击换热性能和平均换热性能计算结果重构成采用最小二乘支持向量机算法建立简单的输入-输出代理模型。

(4)将系统级别数值模拟得到的系统边界条件与小尺度细节代理模型相结合，获得小尺度模型细节在系统边界上的分布。即在空冷单元的迎风面的每个网格上，对网格上的三维风速应用翅片管换热特性代理模型，获得冲击换热系数、平均对流换热系数在迎风面上的二维分布。

采用最小二乘支持向量机将边界条件多变的数值模拟结果得到的翅片管平均对流换热系数hwhole、圆弧段冲击对流换热系数himpine拟合成关于3个来风分速度的函数，从而可将其作为翅片管换热性能的代理模型。

给出训练数据{(xk,yk)}，xk∈Rn，yi∈R，k=1，2，…，M。将训练数据采用非线性函数φ(·):∈Rn→Rnh映射到高维空间，在高维空间寻优进而得到最优回归函数

f(x)=wTφ(x)+b

(14)

该非线性最优回归函数的获得是通过以下二次规划的方法得到的

(15)

式中：c表示对于超过误差范围的样本给予的惩罚因子；ξi为松弛变量；w、b为回归系数向量或系数。

从而得到回归表达式为

(16)

式中K(x,xi)即核函数，为映射函数的点积

K(x,xi)=φ(x)φ(xi)

(17)

本文采用的是RBF径向基核函数，αi、b为拟合得到的系数向量或系数常数。

3 分析及讨论

首先对采用试验关联式的空冷空冷单元不同转速(额定转速69、57.75、34.5、17.25 r/min)下的工况进行了数值模拟，环境温度为5 ℃，从最大转速工况模拟结果中获得距离翅片管束进风口约1 cm的迎风面的三维速度的速度分量的最大变化范围(0～10 m/s)，将各个工况迎风面的三维速度分布提取出来。然后首先对翅片管小尺度模型进行了数值模拟，从数值计算结果中提取出平均对流换热系数和冲击换热系数，并将提取结果用最小二乘支持向量机重构成简单的输入-输出翅片管换热性能代理模型。翅片管的边界条件为：管内壁温度为40 ℃，进口风温为5 ℃，通流风速即z向分速度以及不同的y向、z向分速度，分速度范围在(0～10 m/s)内。最后将不同工况的迎风面三维速度的分布代入翅片管换热性能代理模型中，得到了不同转速下迎风面对应的翅片管束的平均对流换热系数和冲击换热系数分布图。

3.1 最小二乘支持向量机拟合精度

对复杂进风条件下翅片管平均对流传热系数、冲击换热系数的计算结果采用支持向量机进行拟合，表1为拟合精度。由于该算法拟合方法依赖于原始数据，不具有公式一样的广泛推广价值，此处仅给出拟合精度。这种方法能将计算结果很好的程序化，但不具有如试验关联式般明确的物理意义。

表1 支持向量机拟合精度

Table 1 Fitting precision of support vector machine

3.2 复杂进口风条件下的传热系数变化规律

如图3、4为通流风速为2 m/s采用不同的y向、z向分速度的翅片管平均传热系数和圆弧头部最大冲击传热系数的数值模拟结果。可以看出，随着分速度的增加，平均对流换热系数有增大的趋势，但变化幅度不是很大，在图中所示的变化范围内，平均传热系数最大值最小值仅仅相差7 W/(m2K)。而最大冲击换热系数增大的趋势非常显著，最大值最小值相差60 W/(m2K)。

图3 x向通流风速2 m/s时平均换热系数 随y向、z向风速的变化规律

图4 x向通流风速2 m/s时冲击换热系数 随y向、z向风速的变化规律

3.3 空冷单元迎风面风速分布

图5、6分别为空冷单元一侧翅片管束前的迎风面通流风速、绝对风速的二维分布图。对比图5、6，可以看出通流风速、绝对迎面风速在平面上分布很不均匀，且通流风速比绝对迎面风速小很多。通流风速的实验结果可参考文献[13]，与实验结果非常相似，与出口风速相等的通流风速的分布呈现出上半部小下半部大，左右部分相差较大的特征。绝对风速在迎风面上呈现出的上半部小下半部大左右部分相差较大的特征更加明显。分析其原因，风速上半部小下半部大的特征是由于风机出口风受到空冷单元管束的阻挡造成的，左右部分相差较大的特征是带有旋转特性的风机来风造成的。

图5 空冷单元迎风面x方向通流风速分布 Fig.5 Through-velocity distribution at

图6 空冷单元迎风面绝对速度分布

3.4 空冷单元管束迎风面平均对流换热系数二维分布、冲击换热系数二维分布

图7、8为空冷单元一侧翅片管束前面的迎风面对应的翅片管束的平均对流传热系数、冲击换热系数的二维分布图。从这2幅图可以看出，管束迎风面冲击换热系数比平均对流换热系数大很多，最大达到146 W/(m2K)，而平均对流换热系数大部分在30 W/(m2K)左右。这表明空冷单元风机来风对管束的冲击换热远大于管束通流部分的平均对流换热。这是由于管束前端的扁平管圆弧进口段受到了风机来风绝对风速的冲刷，相比通流部分，边界层薄很多，这使得冲击冷却过程的换热系数高出通流部分很多。

图7 空冷单元额定转速管束迎风面平均对流换热系数

图8 额定转速下管束迎风面冲击换热系数分布

3.5 不同转速下冲击换热系数的分布规律

风机来风时，翅片模型圆弧头部处于冲击换热状态，边界层很薄，传热系数很大；迎面风速对管束进口圆弧区域的冲击换热是以风机来风绝对风速进行的；此处风温处接近环境温度，这3个原因导致在冬季工况下，圆弧头部的冲击换热直接影响凝结水的冻结风险的大小。而工程实践中多采用调节风机转速的方法防冻，因此本节对不同风机转速的冲击换热特性进行了数值计算，结果如图9所示。从图9可以知道，随着风机转速的降低，冲击换热系数明显降低。

4 结 论

(1)在空冷风机进风变化范围内，通过翅片管的数值模拟，获得了复杂进风条件下的翅片管平均对流换热系数、圆弧段冲击换热系数，并用支持向量机算法将计算结果简化成输入-输出代理模型。结果表明，进风条件对整体换热能力略有影响，对冲击换热影响较大。

图9 变转速下管束迎风面冲击换热系数分布

(2)通过对空冷单元和翅片空间的多尺度模拟，获得了空冷单元迎风面风速的二维分布，给出了翅片管束的平均对流换热系数、圆弧段冲击换热系数的二维分布。结果表明，平均对流换热系数和冲击换热系数空间分布很不均匀，且冲击传热系数比对流传热系数大一个数量级，增进了对空冷岛防冻机理的理解。

(3)对不同风机转速下的管束进口段的冲击传热性能进行了计算，结果表明，随着转速的下降，冲击传热系数均进一步下降，显示了降低风机转速对降低空冷岛防冻风险的有利影响，有利于下一步同时考虑管内外过程的空冷岛防冻机理的定量分析。

[1]冯丽丽, 杜小泽, 杨勇平, 等.交错蛇形短翅片对扁平管外传热的强化[J].中国电机工程学报, 2012, 32(17): 70-75.Feng Lili, Du Xiaoze, Yang Yongping, et al.Heat transfer enhancement of staggered short wavy fins on flat tube [J].Proceedings of the CSEE, 2012, 32(17): 70-75.

[2]胡汉波.直接空冷式凝汽器翅片散热器流动传热性能及单元流场特性研究[D].重庆大学, 2006.Hu Hanbo.Research on flow and heat transfer of finned heat exchanger of direct-cooled condenser and flow characteristic of condenser unit[D].Chongqing: Chongqing University, 2006.

[3]李书营.蛇形翅片管及空冷单元流动换热规律数值研究[D].北京：北京交通大学, 2009.Li Shuying.Numerical investigation on the heat transfer and fluid flow of wave-finned tube and air-cooled condenser unit[D].Beijing: Beijing Jiaotong University,2009.

[4]侯雪梅.火电厂直接空冷凝汽器蛇形翅片管流动换热性能数值研究[D].济南：山东大学, 2011.Hou Xuemei.Numerical study of finned serpentine tubes’ flow and heat transfer performance in direct air cooling power plant[D].Jinan: Shangdong University, 2011.

[5]马宏涛.空冷凝汽器中流体流动传热特性影响因素的研究[D].北京：北京交通大学,2011.Ma Hongtao.Study on the impact factors of the flow dynamics and heat transfer in the air-cooled condenser[D].Beijing: Beijing Jiaotong University, 2011.

[6]Wen J, Tang D, Wang Z, et al.Large eddy simulation of flow and heat transfer of the flat finned tube in direct air-cooled condensers[J].Applied Thermal Engineering, 2013, 61(2): 75-85.

[7]代善良.迎面风速对直接空冷凝汽器流动换热特性的影响[D].哈尔滨：哈尔滨工业大学,2013.Dai Shanliang.Influence of head wind velocity on flow and heat transfer performance of direct air-cooled condenser[D].Harbin: Harbin Institute of Technology, 2007.

[8]杨立军, 周健, 杜小泽, 等.扁平管外蛇形翅片空间的流动换热性能数值模拟[J].工程热物理学报, 2007, 28(1): 122-124.Yang Lijun, Zhou Jian, Du Xiaoze, et al.Numerical analysis of flow and heat transfer outside the flat wave-finned tube[J].Journal of Engineering Thermophysics, 2007, 28(1): 122-124.

[9]张凯峰, 杨立军, 杜小泽, 等.空冷凝汽器波形翅片扁平管管束外空气流动传热特性[J].中国电机工程学报, 2008, 28(26): 24-28.Zhang Kaifeng, Yang Lijun, Du Xiaoze, et al.Flow and heat transfer characteristics of cooling air outside the flat wave-finned tube banks of air-cooled condensers[J].Proceedings of the CSEE, 2008, 26: 006.

[10]冯丽丽, 杜小泽, 杨勇平, 等.扁平管交错蛇形短翅片强化传热的研究[J].工程热物理学报, 2012, 33(7): 1201-1204.Feng Lili, Du Xiaoze, Yang Lijun, et al.Study on flow and heat transfer of staggered short wavy fins on flat tube[J].Journal of Engineering Thermophysics, 2012, 33(7): 1201-1204.

[11]Wen Juan, Tang Dawei.Large eddy simulation of flow and heat transfer of the flat finned tube in direct air-cooled condensers[J].Applied Thermal Engineering, 2013, 61(2): 75-85.

[12]Zhang Hui, Yue Yujin, Du Xiaoze, et al.Numerical investigation of flow fields of air-cooled condenser cell including fan blade passages[J].HEFAT 2011, 2011.

[13]卜永东.直接空冷单元流动传热特性及空气流场优化组织[D].北京：华北电力大学,2013.Bo Yongdong.Flow and heat transfer characteristic and air-side flow field optimization of air-cooled condenser cell[D].Beijing: North China Electric Power University, 2013.

[14]曾时明.电站直接空冷式翅片管换热和流动规律及空冷单元流场特性数值研究[D].北京: 北京交通大学, 2007.Zeng Shiming.Numerical investigation of the heat transfer and fluid flow of finned-tube and the flow characteristics of a condenser unit in direct air-cooled power plant[D].Beijing: Beijing Jiaotong University, 2007.

[15]杨立军, 杜小泽, 张辉, 等.电站空冷凝汽器轴流风机阵列集群效应的数值研究[J].科学通报, 2011, 56(15): 1232-1239.Yang Lijun, Du Xiaoze, Zhang Hui, et al.Numerical investigation on the cluster of an array of axial flow fans for air-cooled condensers in a power plant[J].Chinese Sci Bull, 2011, 56(15): 1232-1239.

[16]Li Li, Du Xiaoze, Yang Lijun, et al.Numerical simulation on flow and heat transfer of fin structure in air-cooled heat exchanger[J].Applied Thermal Engineering, 2013, 59(1): 77-86.

[17]杜小泽,杨立军, 金衍胜, 等.火电站直接空冷凝汽器传热系数实验关联式[J].中国电机工程学报，2008,14(28)：32-37 Du Xiaoze, Yang Lijun, Jin Yansheng, et al.Development of experimental correlation for heat transfer coefficient of direct air-cooling condenser in power plant[J].Proceedings of the CSEE, 2008, 14(28):32-37

[18]肖烈晖, 杜小泽, 杨立军，等 .环境风作用下2×1 000 MW直接空冷机组空冷岛布局[J].电力建设, 2015, 36(6): 7-13.Xiao Liehui, Du Xiaoze, Yang Lijun，et al .Air-cooling condenser layout of 2×1 000 MW direct air-cooled units under effect of ambient natural wind[J].Electric Power Construction, 2015, 36(6): 7-13.

[19]魏高升, 陈林, 李兵，等 .板式蒸发式冷凝器热质传递过程模型分析[J].电力建设, 2014, 35(6): 18-21.Wei Gaosheng, Chen Lin, Li Bing，et al.Model analysis on heat and mass transfer process in plate-type evaporative condenser[J].Electric Power Construction, 2014, 35(6): 18-21.

(编辑:蒋毅恒)

CHENG Tongrui, DU Xiaoze, YANG Lijun

(Key Laboratory of Condition Monitoring and Control for Power Plant Equipment of

Ministry of Education, North China Electric Power University, Beijing 102206, China)

For the antifreezing industrial demand of direct air-cooling condenser in winter, the model on finned tube scale was coupled with that on air-cooling unit scale, both the whole heat transfer characteristics in finned tube bundle and the impingement heat transfer characteristic of condensate water in air-cooling finned tube were studied, which could provide the third thermal boundary condition for the calculation of the in-tube condensation heat transfer and freezing in direct air-cooling condenser under the condition of low temperature operation in winter.The air-cooling condenser cell model and finned tube model were constructed.Then, the simulation results of finned tube under 3D inlet wind velocity boundary condition was simplified as input-output type agent model.The velocity distribution at the windward surface and symmetric boundary was obtained from air-cooling condenser cell numerical simulation.Then the finned tube agent model was applied in the symmetric boundary to obtain the spatial distribution of the impingement heat transfer characteristics and average convective heat transfer characteristics of finned tube at symmetric boundary.The numerical calculation results show that, the impinging heat transfer coefficient has a larger magnitude than the average convective heat transfer coefficient; at the windward surface, the distributions of those two coefficients are both up-down asymmetrical and left-right asymmetrical; and the impinging heat transfer coefficient decreases with the decrease of the rotational speed.

direct air-cooling condenser; cross-scale simulation; heat transfer; freezing

国家重点基础研究发展计划项目(973项目)(2015CB251503)。

TM 621; TK 124

A

1000-7229(2015)08-0015-07

10.3969/j.issn.1000-7229.2015.08.003

2015-05-25

2015-07-10

程通锐(1984)，男，博士研究生 主要研究方向为电站优化与强化传热；

杜小泽(1970)，男，通信作者，博士生导师，主要研究方向为电站优化与强化传热；

杨立军(1971)，男，博士生导师，主要研究方向为电站优化与强化传热。

Project Supported by National Basic Research Program of China (973 Program)(2015CB251503).