蒸汽冷凝器流动传热过程的准三维CFD模拟算法
吴家豪1, 代守宝2, 张归华1, 王雄师1, 赵艳伟1, 吴玉新1    
1. 清华大学 能源与动力工程系, 热科学与动力工程教育部重点实验室, 北京 100084;
2. 中国船舶集团有限公司第七〇三研究所, 哈尔滨 150078
摘要:管侧冷却水在流向和横向(垂直于流向)上的分布对管壳式蒸汽冷凝器内部各物理量的分布和冷凝器的设计优化都会产生影响, 因此有必要对其进行研究。该文提出了一种考虑冷却水横向分布的准三维CFD模拟算法, 并基于自主开发的二维CFD冷凝过程数值计算代码ConDesign-2D对已有算例进行了模拟。结果表明, 冷却水沿流向的温升不容忽略, 它直接降低了换热温差, 进而影响到了冷凝速率等物理量的分布, 冷却水的横向分布也使得冷凝速率分布更加均匀。相比二维模拟, 该算法能解析出流场在冷却水流向上的第三维分布信息, 又避免了全三维模拟带来的高昂计算成本, 可以应用于实际冷凝器的模拟计算。
关键词准三维    冷却水分布    冷凝器    计算流体动力学(CFD)模拟    
Quasi-3D CFD algorithm for the flow and heat transfer process in steam condensers
WU Jiahao1, DAI Shoubao2, ZHANG Guihua1, WANG Xiongshi1, ZHAO Yanwei1, WU Yuxin1    
1. Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Department of Energy and Power Engineering, Tsinghua University, Beijing 100084, China;
2. No. 703 Research Institute of China State Shipbuilding Co., Ltd., Harbin 150078, China
Abstract: [Objective] Shell-and-tube steam condensers are widely used in power stations and ships; their performance, which requires an in-depth understanding of the internal flow and heat transfer process, ensures device efficiency and safety. The distribution of tube-side cooling water in its flow (z-direction) and transverse (x-y plane) directions influences the distributions of various physical quantities inside a condenser and its design optimization, which cannot be calculated by purely two-dimensional (2D) computational fluid dynamics (CFD) simulation. Moreover, the complexity of the two-phase flow, turbulence, phase change, and heat transfer mechanisms inside the condenser makes full 3D CFD simulations computationally expensive. [Methods] Herein, a quasi-3D algorithm considering the transverse distribution of cooling water was proposed and applied to simulate the condensation process in a shell-and-tube steam condenser. The tube bundle region was simplified through porous media assumption, introducing extra resistance source terms to momentum transport equations. The condensation source term in the continuity equation was computed based on equivalent thermal resistance, referring to the summation of cooling-water convection, tube wall, condensate water, and noncondensable gas (air) thermal resistance. Considering the temperature rise along the z-direction, the quasi-3D algorithm split the condenser into several sections along the z-direction, performed 2D simulation in the midplane of each section, and used the secant iteration method to balance the steam pressure drop. The simulation was conducted based on ConDesign-2D, a self-developed 2D CFD code for the condensation process of shell-and-tube condensers, which adoptedunstructured meshes and the collocated-grid-based SIMPLE algorithm. [Results] The temperature rise of cooling water along the flow direction directly reduced the heat transfer temperature difference by 10% and affected the distribution of important physical quantities such as the condensation rate. Compared with the reference condensation rates of different sections, the calculated values considering the transverse distribution of cooling water are more accurate than those that do not. Additionally, the transverse distribution of the cooling water led to a more uniform condensation rate distribution due to the negative feedback between cooling-water temperature and condensation rate. However, this variable had less influence on the distribution of noncondensable gas (air). The comparison with 2D results revealed that 2D and quasi-3D simulations gave similar results on the midplane of the condenser, which illustrates the linearity of field distributions in the z-direction. According to the computation complexity analysis, the complexity of the full 3D simulation can be 2-3 orders of magnitude higher than that of 2D simulation, whereas quasi-3D simulation can be only 1 order of magnitude higher. [Conclusions] Compared with 2D simulation, the proposed quasi-3D algorithm can compute the 3D distribution information of the flow field in the cooling-water flow direction and avoid the high computational cost brought by full 3D simulation. Therefore, for the rapid and accurate prediction of 3D flow and the heat transfer process inside a practical shell-and-tube condenser, the proposed quasi-3D algorithm is preferred over 2D and full 3D simulations.
Key words: quasi-3D    cooling-water distribution    condenser    computational fluid dynamics (CFD) simulation    

管壳式蒸汽冷凝器被广泛用于电站和船舶等,其性能对于保障装置热效率及安全性具有重要意义,为此需要深入了解其内部流动和传热过程,以便对其做进一步的优化设计。计算流体动力学(computational fluid dynamics,CFD)模拟是研究冷凝器内部流动和传热过程的重要辅助手段,并日益被广泛应用。然而,冷凝器内部存在三维流动、两相流、湍流、相变和传热等多个复杂的物理机制,使得对其进行完全的、三维的CFD模拟的计算成本高昂。为此,迫切需要在冷凝器设计过程中建立简化的CFD模拟方法,从而能够在符合实际物理过程的基础上,有效降低计算成本并准确预测冷凝器性能及传热过程。

在冷凝器数值模拟的模型简化方面,为降低管束区的计算代价,Patankar等[1]提出分布阻力方法,忽略冷却管的固体壁面,将管束区视为多孔介质,并用孔隙率来衡量冷却管的密集程度,用阻力源项和冷凝源项来衡量管束带来的流动阻力和管束区的气体冷凝速率。目前,国内外大多数冷凝器数值模拟研究都会使用多孔介质假设以降低计算成本。

冷凝器数值模拟的空间维度主要分为一维、二维、准三维和全三维。其中,一维模拟只考虑参数在单一方向上的变化,准确性低,经验性强,因此应用范围相对有限,主要用于估算整体传热特性[2]。二维模拟的应用较多[3-6],只在垂直于冷却水流向的二维平面上进行模拟,其合理性在于冷凝器壳侧的流动主要都沿着此平面方向。准三维模拟是在考虑冷却水温升的基础上,在冷却水流向不同位置处进行二维模拟,可以进一步解析出冷凝器内各物理量的三维分布信息[7-8]。全三维模拟的出现源于计算机性能的提升[9],它最符合物理实际,但计算成本较高,且在管束和壳体优化设计的指导效果上并无明显优势。因此,二维和准三维模拟仍然是冷凝器壳侧流场模拟的首选。最后需要指出的是,冷凝器常使用隔板来支撑管束,这会限制壳侧沿第三个方向的流动,这进一步说明了二维模拟的合理性;隔板将冷凝器沿冷却水流向分为若干个扇区,且隔板之间的间距相对较小,这进一步说明了准三维模拟的合理性,使得工程师可以将每个扇区视为没有厚度的平面来进行二维模拟,从而得到准三维模拟的结果。

准三维模拟中的冷却水温升是指冷却水温度在流向上的变化或分布,而冷却水在横向(垂直于流向)上的分布对于计算结果和冷凝器设计优化的影响却未得到足够的重视。冷却水的横向分布主要包括2点:流量的分布和温度的分布。前者主要由冷却管的排布决定,管路越密集的地方,冷却水流量就越高,其冷却效果也越好;而后者则直接影响对前者的优化设计,例如冷却水温度较高的地方说明其换热和冷凝较强,因此可以考虑在对应的地方布置更为密集的冷却管。另外,冷却水的横向分布会影响到其他物理量的分布,例如换热温差的分布改变将直接影响冷凝速率的分布情况,进而影响蒸汽凝结量和不凝性气体的分布情况,这都会影响冷凝器的设计优化。因此,有必要对冷却水的横向分布进行研究。

本文提出了一种考虑冷却水横向分布的准三维CFD模拟算法,在简单二维模拟的基础上,考虑了冷却水沿流向的温升以及流向不同位置处的蒸汽压降平衡,并基于割线法对蒸汽入口速度分布进行迭代求解,使算法既能模拟出冷凝器内各物理量的三维分布信息,又能模拟出冷却水的横向分布及其对其他物理量分布所产生的影响,且避免了全三维模拟的高昂计算成本,兼顾了计算精度和效率,有利于提升实际冷凝器的设计优化效果。

1 模型与算法 1.1 基本假设

本文在考虑冷凝器壳侧流动特点的基础上,进行了以下简化假设:1)将管束区等效为多孔介质,使用多孔介质模型进行计算;2) 忽略冷凝液占据的体积及其与混合气体间的相互作用,只计算混合气体的单相流动,因此忽略气液速度差造成的动量方程源项;3) 将组分复杂的空气简化为纯净气体,即壳侧只有空气和蒸汽2种组分;4) 认为冷凝液和蒸汽始终处于饱和状态,因此壳侧气体的温度可由压力查表得来;5) 将冷凝器沿冷却水流向划分为若干扇区,将每个扇区视为没有厚度的平面来进行二维模拟,并考虑扇区之间的相互作用,由此进行准三维模拟。

1.2 控制方程

质量方程的Einstein求和形式为

$ \frac{\partial}{\partial x_j}\left(\varepsilon \rho u_j\right)=-S_{\mathrm{m}} . $ (1)

其中:j为求和哑指标,j=1, 2;x1x2分别表示空间坐标xyu1u2为对应的速度分量;ε为孔隙率;ρ为混合气体的密度;Sm为冷凝源项,表示单位体积冷凝速率,单位为kg·m-3·s-1。动量方程的求和形式为

$ \begin{gathered} \frac{\partial}{\partial x_j}\left(\varepsilon \rho u_j u_i\right)-\frac{\partial}{\partial x_j}\left(\varepsilon \mu_{\mathrm{eff}}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\right)= \\ -\varepsilon \frac{\partial p}{\partial x_i}-S_{\mathrm{m}} u_i+S_{u, i}, i=1, 2 . \end{gathered} $ (2)

其中:i为1和2分别表示xy维度,p为压力,μeff为等效黏度,Su, i为管束排布造成的阻力源项。

组分方程只包括空气的质量分数输运方程:

$ \frac{\partial}{\partial x_j}\left(\varepsilon \rho u_j Y_{\text {air }}\right)=\frac{\partial}{\partial x_j}\left(\varepsilon D_{\text {eff }} \frac{\partial Y_{\text {air }}}{\partial x_j}\right) . $ (3)

其中:Yair为空气质量分数,Deff为等效扩散系数。

最后,补充性的方程包括蒸汽的质量分数计算式以及基于理想气体混合律的混合气体密度计算式:

$ Y_{\mathrm{s}}=1-Y_{\text {air }} \text {, } $ (4)
$ \frac{1}{\rho}=\frac{Y_{\text {air }}}{\rho_{\text {air }}}+\frac{Y_{\mathrm{s}}}{\rho_{\mathrm{s}}} \text {. } $ (5)

其中:Ys为蒸汽的质量分数,ρairρs分别为空气和蒸汽的密度。

1.3 阻力源项

Su, i[4]

$ S_{u, i}=\varepsilon \xi_i \rho u_i|\boldsymbol{U}| . $ (6)

其中:ξi为压力损失系数,U为速度矢量。

$ \xi_i=2 \frac{f_i}{P_{\mathrm{t}}}\left(\frac{P_{\mathrm{t}} \varepsilon}{P_{\mathrm{t}}-D_{\mathrm{o}}}\right)^2 ; $ (7)
$ f_i= \begin{cases}0.619 Re_i^{-0.198}, & Re_i \leqslant 8000 ; \\ 1.156 Re_i^{-0.2647}, & 8000<Re_i \leqslant 2 \times 10^5 ;\end{cases} $ (8)
$ R e_i=\frac{\rho u_i D_0}{\mu} . $ (9)

其中:Pt为管心距,Do为管外径,Rei为速度ui的Reynolds数。

1.4 冷凝源项

壳侧蒸汽冷凝释放的热量全部用于加热冷却水,由此可得

$ S_{\mathrm{m}}=\frac{A \Delta T}{R_{\text {total }} h_{\mathrm{lg}} V} . $ (10)

其中:AV分别为控制体中的管束表面积和控制体体积,ΔT=Ts-Tcw为蒸汽与冷却水的温差,hlg为蒸汽的冷凝潜热,Rtotal为蒸汽与冷却水之间的总热阻。Rtotal由冷却水对流热阻Rcw、管壁导热热阻Rtb、液膜热阻Rc和不凝气气膜热阻Ra构成:

$ R_{\text {total }}=\frac{D_{\mathrm{o}}}{D_{\mathrm{i}}} R_{\mathrm{cw}}+R_{\mathrm{tb}}+R_{\mathrm{c}}+R_{\mathrm{a}} . $ (11)

其中Di为管内径。

Rcw使用Petukhov关联式[10]表示:

$ \begin{aligned} & \frac{1}{R_{\mathrm{cw}}}=\frac{\lambda_{\mathrm{cw}}}{D_{\mathrm{i}}} Re_{\mathrm{cw}} Pr_{\mathrm{cw}}\left(f_R / 8\right) /\left[1+3.4 f_R+\right. \\ & \left.\left(11.7+1.8 Pr_{\mathrm{cw}}^{-1 / 3}\right)\left(f_R / 8\right)^{1 / 2}\left(Pr_{\mathrm{cw}}^{2 / 3}-1\right)\right], \end{aligned} $ (12)
$ f_R=\left(0.79 \ln R e_{\mathrm{cw}}-1.64\right)^{-2} . $ (13)

其中:λcwRecwPrcw分别为冷却水的热导率、Reynolds数和Prandtl数。

$ R_{\mathrm{tb}}=\frac{D_{\mathrm{o}} \ln \left(D_{\mathrm{o}} / D_{\mathrm{i}}\right)}{2 \lambda_{\mathrm{tb}}} . $ (14)

其中λtb为管壁热导率。

Rc使用Hu等[11]的管高修正模型表示:

$ R_{\mathrm{c}}=R_{\mathrm{c}}^{\prime} n_{\mathrm{tb}}^{0.33 Y} . $ (15)

其中:ntb为从高往低数的冷却管排数;Yy的线性函数,在最高排管处为0,最低排管处为1;R′c为单管液膜热阻。

$ R_{\mathrm{c}}^{\prime}=\left[N u_{\mathrm{c}} \frac{\lambda_{\mathrm{c}}}{D_{\mathrm{o}}}\left(1+0.0095 R e_{\mathrm{s}}^{11.8 / \sqrt{N u_{\mathrm{c}}}}\right)\right]^{-1} . $ (16)

其中:λc为冷凝液的热导率,Res为蒸汽的Reynolds数,Nuc为液膜Nusselt数。

$ N u_{\mathrm{c}}=0.725\left[\frac{D_{\mathrm{o}}^3 g \rho_{\mathrm{c}}^2 h_{\mathrm{lg}}}{\mu_{\mathrm{c}} \lambda_{\mathrm{c}}\left(T_{\mathrm{s}}-T_{\mathrm{cw}}\right)}\right]^{\frac{1}{4}} . $ (17)

其中ρcμc分别为冷凝液的密度和黏度。

Ra使用Berman等[12]给出的经验关联式表示:

$ \begin{gathered} \frac{1}{R_{\mathrm{a}}}=a \frac{D_{\mathrm{k}}}{D_{\mathrm{o}}} Re_{\mathrm{s}}^{1 / 2}\left(\frac{p}{p-p_{\mathrm{s}}}\right)^b p^{1 / 3} \cdot \\ \left(\frac{\rho_{\mathrm{s}} h_{\mathrm{lg}}}{T_{\mathrm{s}}}\right) \frac{1}{\left(T_{\mathrm{s}}-T_{\mathrm{cw}}\right)^{1 / 3}} . \end{gathered} $ (18)

其中:Dk为扩散系数,ps为蒸汽分压,ab为系数。

$ (a, b)= \begin{cases}(0.52, 0.7), & Re_{\mathrm{s}} \leqslant 350 ; \\ (0.82, 0.6), & Re_{\mathrm{s}}>350 .\end{cases} $ (19)
1.5 准三维算法

将冷凝器沿冷却水流向平均分为N个扇区,在每个扇区的中截面上进行二维CFD模拟,从而得到冷凝器各参数在冷却水流向这一维度上的大致变化情况。在N个扇区的耦合方面,准三维算法主要考虑了2个因素:冷却水沿流向的温升和N个扇区的蒸汽压降平衡。

1.5.1 冷却水温升

对每个扇区截面,认为蒸汽冷凝释放的热量全部用于加热冷却水,因此有

$ \Delta T_{\mathrm{cw}}=\frac{S_{\mathrm{m}}(L / N) A_{\text {cross }} h_{\mathrm{lg}}}{M_{\mathrm{cwc}} c_{\mathrm{p}, \mathrm{cw}}} . $ (20)

其中:ΔTcw为所考虑区域内冷却水沿流向的温升,L为冷却管总长度,cp, cw为冷却水比热容,SmAcrossMcw分别为所考虑区域的平均冷凝源项、横截面积和冷却水质量流量。所考虑区域在模拟时有2种配置:单个的计算网格,在管束区内包含所有计算网格的全计算域。其中前者可以计算出冷却水的横向分布,包括温度和质量流量分布;而后者只能计算出冷却水物理量在横向平面的平均值。

对所考虑区域,认为冷却水温度Tcw(取中截面处的冷却水温度)和冷却水出口温度Tcw, out

$ \left\{\begin{array}{l} T_{\mathrm{cw}, n}=T_{\mathrm{cw} \cdot \text { in }, n}+\Delta T_{\mathrm{cw} \cdot n} / 2, \\ T_{\mathrm{cw}, \text { out }, n}=T_{\mathrm{cw} \cdot \text { in }, n}+\Delta T_{\mathrm{cw}, n}, n=1, 2, \cdots N . \end{array}\right. $ (21)

其中n为扇区编号。冷却水入口温度Tcw, in,1为给定值,Tcw, in, n (n=2, 3, …, N)等于上一扇区的冷却水出口温度,即

$ T_{\mathrm{cw}, \text { in }, n}=T_{\mathrm{cw}, \text { out }, n-1}, n=2, 3, \cdots N \text {. } $ (22)

因此,冷却水温升算法是一种顺序计算算法,需要在扇区n的二维CFD模拟完成后才能对扇区(n+1)进行计算。

1.5.2 蒸汽压降平衡

如果N个扇区的壳侧入口流量相同,那么在冷却水温升的影响下,各扇区的冷凝量会有所差异,进而造成各扇区的沿程压降有所不同。然而,在实际冷凝器运行过程中,由于N个扇区的壳侧入口和出口都分别是相通的,因此它们的压降应该相等。因此,需要设计算法,通过调整N个扇区的入口流量,来使它们的压降相平衡,得到符合物理实际的准三维模拟结果。

由于本文计算时采用压力出口边界条件,N个截面的出口压力相同,因此算法目标即为N个截面的入口压力相等。考虑冷却水温升和蒸汽压降平衡的准三维模拟算法流程如图 1所示。其中:maxiter为最大迭代次数,tol为压差收敛判据值,φ为迭代的松弛因子,k为当前迭代数,vn (n=1, 2, …, N)为不同扇区的入口蒸汽速度,v0为给定平均速度值,α为入口速度分布斜率参数。

图 1 准三维模拟算法流程

算法进行了3项假设:1)认为壳侧入口流量沿z方向(冷却水流向)线性分布,并用参数α1来间接衡量其分布的斜率;2)认为入口压力沿z方向单调分布,进而可以将扇区1和N的入口压力之差Δp=p1-pN作为最小化目标;3)认为Δpα1之间的线性关系较好,因此可以使用割线法来对此优化问题进行迭代求解。另外,当N取1时即为二维模拟,此时无需进行蒸汽压降平衡的迭代。

2 研究对象与计算程序 2.1 研究对象

本文模拟的算例为文[13]中的实验室级冷凝器,具体结构如图 2所示,管束区内有20×20根的叉排冷凝管束,冷却水流向垂直于蒸汽流动的方向。

图 2 冷凝器结构与几何尺寸[13]

算例的其他几何参数与边界条件如表 1所示。

表 1 算例几何参数与边界条件[13]
类别 参数 数值
几何参数 Do/mm 25.4
管壁厚δ/mm 1.25
Pt/mm 34.9
边界条件 蒸汽入口流量Min/(kg·s-1) 2.032
蒸汽出口压力pout/Pa 27 530
冷却水入口流速ucw, in/(m·s-1) 1.19
Tcw, in/℃ 17.8

表 1可以计算得到本算例的孔隙率为

$ \varepsilon=\left\{\begin{array}{l} 1, \\ 1-\frac{\pi}{2 \sqrt{3}}\left(\frac{D_{\mathrm{o}}}{P_{\mathrm{t}}}\right)^2=0.5196 . \end{array}\right. $ (23)
2.2 计算程序

本文使用自主开发的二维CFD冷凝过程数值计算代码ConDesign-2D对所选算例进行模拟计算。代码的核心求解算法基于有限体积法,网格使用非结构网格,以适应实际的集成冷凝器的几何复杂性。梯度计算基于Gauss-Green方法,对非结构网格的非正交部分,采用显式延迟修正的方式处理[14]。在压力—速度耦合方面,使用基于同位网格的SIMPLE算法以解决非结构网格无法使用交错网格的问题。基于开源程序bamg进行网格的生成。最终使用的计算网格如图 3所示,网格单元数为1 824。

图 3 计算网格示意图

3 结果与讨论

对于冷却水在流向上的参数变化,将二维模拟和准三维模拟的结果进行对比;对于冷却水在横向扇区截面上的参数变化,如1.5节所述,将冷却水温升算法中两种不同所选区域下的计算结果进行对比,即使用式(20)计算温升时,平均冷凝源项、横截面积和冷却水质量流量的考虑区域分别选择最小的控制体和含冷却水管的全计算域。本文准三维模拟的扇区截面数设为5。部分计算结果与文[13]进行了对比,以验证本文计算的准确性。

考虑冷却水横向分布时,z方向不同位置处的冷却水温度分布情况如图 4所示。可以看出,在z方向上,冷却水温度不断升高,到出口处的平均升温为5.26℃,这会直接降低壳侧蒸汽与冷却水之间的温差,进而降低冷凝率。定量地,壳侧蒸汽温度约为67.3℃,入口处的换热温差为67.3℃-17.8℃=49.5℃,因此出口处的平均换热温差降低了约5.26/49.5=10.6%。这一数值直接说明了不应忽略冷却水温度在流向上的变化,也初步说明了准三维模拟的必要性。在x-y平面上,冷却水温度分布不均匀:在靠近出口处较低,而在靠近管束区边缘处较高,这是因为冷却水温度与冷凝速率直接相关,蒸汽刚进入管束区时的冷凝速率高,换热量大,因此冷却水升温高,在出口处则反之。定量地,冷却水到出口处的最低升温和最高升温分别为0.92和8.79℃,对应的平均换热温差分别为48.58和40.71℃,和出口处的平均换热温差44.24℃分别相差9.8%和-9.2%。这一数值说明不应忽略冷却水的横向温度分布,特别是在冷却水下游处,因为由图 4可以看出冷却水温度的不均匀性沿流向逐渐增大。另外,本文算法和代码中也考虑了冷却水流量的横向分布,这对于冷却管非均匀排布的冷凝器计算十分必要,但是由于本文所选算例的冷却管排布均匀,因此未予以研究。

图 4 Tcw在不同z处的分布

蒸汽压降平衡算法迭代时,所有的计算都在3个迭代步后收敛,这是割线法的最小迭代次数,说明前述“Δpα1之间的线性关系较好”的假设成立。最终收敛时,不同扇区截面的冷凝率如图 5所示。可以看出,考虑冷却水横向分布的结果比不考虑的更接近文献值,两者与文献值的最大相对偏差分别为0.16%和0.75%,进一步说明考虑冷却水横向分布的计算结果更为准确。

图 5 不同扇区的冷凝率

图 6对比了二维模拟、准三维模拟和文[13]中的不同高度处的热流密度变化情况。热流密度越大,表示换热强度越强。z方向不同位置处的各热流密度曲线有所差异,这进一步说明了三维效应对冷凝器内换热过程的影响:对于第3排管,各x处的换热强度都沿z方向逐渐增强,而对于第8、13和18排管,则是靠近管束区边缘处的换热强度沿z方向逐渐增强,而位于管束区中间位置的则逐渐减弱。定量地,第3、8、13和18排管在z为0.9L处的热流密度相比z为0.1L处的最大相对偏差分别为13.6%、14.1%、14.5%和14.8%。另外,二维模拟和准三维模拟在z为0.5L处的结果几乎相同,这说明大部分参数沿z方向分布的线性较好,使得根据细分扇区进行顺序计算的准三维模拟结果与将所有扇区进行整体考虑的二维模拟结果相差不大。需要指出的是,尽管如此,二维模拟却只能得到z为0.5L处的结果,不能像准三维模拟一样能解析出流场在冷却水流向上的第三维分布信息。最后,图中模拟值与实验值的差异主要是由于模型的各简化假设造成的,例如将两相流简化为单相流等,这说明本文所用模型有待进一步改善。

图 6 不同高度处的热流密度分布

图 7进一步定性地展示了不同模拟设置下和z方向不同位置处的冷凝速率分布。冷凝速率越大,表示换热强度越强。首先,如前所述,各情况下冷凝速率的分布特点相同:蒸汽刚进入管束区时的冷凝速率高,而在靠近出口处,由于水蒸气量减少、Reynolds数降低和总热阻增大等多方面原因,导致冷凝速率降低。其次,对比图 7a7c可知,有两点结论与前述图 6结论一致,即关于不同高度处换热强度沿z方向的变化情况的结论,以及二维模拟与准三维模拟在z为0.5L处的结果差异很小的结论。另外由云图可以进一步发现,沿z方向冷凝速率或热流密度分布的均匀性是逐渐增大的。最后,对比图 7b7c可知,考虑冷却水横向分布会使得计算得到的冷凝速率分布更均匀,尤其是在下游,这是因为它考虑了冷凝速率和冷却水温度之间的耦合作用:冷凝速率高的地方换热强,会使得冷却水升温高,而这会降低换热温差,从而降低冷凝速率,所以最终表现出来的结果是冷凝速率在其较高的地方会有所降低,而在较低的地方则会有所升高,使得最终分布比不考虑冷却水横向分布的结果更均匀。定量地,不考虑和考虑冷却水横向分布在z为0.9L处冷凝速率的最大值分别为5.95和5.56 kg/(m3·s),后者比前者小6.6%;最小值分别为0.76和0.82 kg/(m3·s),后者比前者大7.9%。

图 7 不同模拟设置所计算得到的Sm分布

由于在冷凝器的工业设计过程中,不凝结空气聚集点的预测是最重要的环节之一,它有利于管束分布的设计优化,因此将不同计算设置下和z方向不同位置处的空气质量分数分布展示于图 8。首先,空气质量分数沿z方向逐渐降低,这是因为沿z方向换热温差降低,导致冷凝速率和蒸汽凝结量降低,所以空气占比增大。其次,二维模拟与准三维模拟在z为0.5L处的空气分布仍然差异很小,这与前述结论一致。最后,考虑冷却水横向分布对空气质量分数分布的影响不大。定量地,不考虑和考虑冷却水横向分布在z为0.9L处的空气质量分数最大值分别为0.0168和0.0172,后者比前者高2.38%。

图 8 不同模拟设置所计算得到的Yair分布

最后,表 2对比了二维模拟、准三维模拟和全三维模拟的计算复杂度、计算时间等信息。可以看出,虽然全三维模拟具有很高的精度,能解析三维信息,但其计算复杂度较高,一般会比二维模拟高2到3个数量级,因此无法快速给出冷凝器的流场结果。另一方面,二维模拟虽然计算较快,但无法得到流场的三维分布。相比之下,准三维模拟不仅可以解析三维流场信息,而且其计算复杂度和计算时间只是二维模拟的常数倍,一般只高出1个数量级。综合来看,在冷凝器的三维流动传热过程快速计算方面,准三维模拟要优于二维和全三维模拟。

表 2 二维、准三维和全三维模拟的对比
项目 二维模拟 准三维模拟 全三维模拟
计算复杂度 O(m2) O(m2NK) O(m3)
计算时间/s 200 200NK >>200NK
计算精度 中截面处和准三维相同 较高 最高
第三维信息
注:m为平均到每个空间维度上的网格数;K为准三维模拟实际迭代数;NK在本文中分别为5和3。

4 结论

本文建立了蒸汽冷凝器流动传热过程的准三维CFD模拟算法,在算法中考虑了冷却水的流向温升和横向分布,并使用自主开发的基于非结构网格和有限体积法的二维CFD冷凝过程数值计算代码ConDesign-2D对参考算例进行模拟,可以在较低的计算成本下,计算出冷凝器内部流动和传热过程的三维信息、冷却水的横向分布及其对其他物理量分布所产生的影响。

二维模拟与准三维模拟的结果,以及考虑和不考虑冷却水横向分布的结果对比表明:首先,冷却水沿流向的温升不容忽略,它直接降低了换热温差,进而影响到了热流密度、冷凝速率和空气质量分数等物理量沿冷却水流向上的分布情况;其次,冷却水的横向分布也应当考虑,此时的计算结果更接近文献值,而且它影响到了各物理量的横向分布情况,例如使得冷凝速率分布更加均匀;最后,准三维模拟与二维模拟在中截面上的模拟结果差异较小,但是后者不能像前者一样解析出流场在冷却水流向上的第三维分布信息,而且这一结论对于各参数沿冷却水流向非线性分布的冷凝器不一定成立。

为进一步研究冷却水对冷凝器内部流动传热过程的影响,未来的工作方向包括:1)使用更准确的冷却水物性。目前的计算使用的是冷却水在其平均温度下的物性参数,而实际上其物性是随着其温度变化的,对于本算例的冷却水最大温差8.79℃以及其他冷凝器下更大的温差,其物性的变化可能会对结果造成一定程度上的影响。2) 对其他实际冷凝器算例进行研究,尤其是对于冷却管排布不均匀的冷凝器,研究是否考虑冷却水流量横向分布对结果的影响。3) 对本算例的冷却管排布进行优化,例如可以在不改变总管数的前提下,让管束区边缘的冷却管更密集,计算此时冷凝器在换热性能、凝结比等方面的提升程度。

参考文献
[1]
PATANKAR S V, SPALDING D B. A calculation procedure for the transient and steady-state behaviour of shell-and-tube heat exchangers[M]. London: Imperial College of Science and Technology, Department of Mechanical Engineering, 1972.
[2]
BECKETT G, DAVIDSON B J, FERRISON J A. The use of computer programs to improve condenser performance[C]//Condensers: Theory and Practice. Symposium. Oxford, UK: Pergamon Press, 1983: 97-110.
[3]
俞茂铮, 姚秀平, 汪国山, 等. 大功率汽轮机凝汽器汽相流动与传热特性的数值分析[J]. 动力工程, 1995, 15(6): 42-49.
YU M Z, YAO X P, WANG G S, et al. Numerical analysis of the steam flow behavior and heat transfer performance of condensers of large steam turbines[J]. Power Engineering, 1995, 15(6): 42-49. (in Chinese)
[4]
RHODES D B, CARLUCCI L N. Predicted and measured velocity distributions in a model heat exchanger[C]//Proceedings of the International Conference on Numerical Methods in Nuclear Engineering. Montreal, Quebec: American Nuclear Society/Canadian Nuclear Society, 1983: 935-948.
[5]
NUMRICH R, BERKEMEIER D, RENNHACK R. Heat and mass-transfer on condensation of water-vapor in the presence of air at pressure up to 21-bar[J]. Chemie Ingenieur Technik, 1989, 61(8): 666-667. DOI:10.1002/cite.330610833
[6]
ZHANG C, SOUSA A C M, VENART J E S. The numerical and experimental study of a power plant condenser[J]. Journal of Heat Transfer, 1993, 115(2): 435-445. DOI:10.1115/1.2910696
[7]
姚秀平, 李红梅, 朱光宇, 等. 电站凝汽器汽相流动与传热特性的准三维数值分析[J]. 西安交通大学学报, 1998, 32(12): 29-33.
YAO X P, LI H M, ZHU G Y, et al. Quasi-three-dimensional analysis of steam flow and heat transfer of a large power plant condenser[J]. Journal of Xi'an Jiaotong University, 1998, 32(12): 29-33. (in Chinese)
[8]
ZHANG C, ZHANG Y. A quasi-three-dimensional approach to predict the performance of steam surface condensers[J]. Journal of Energy Resources Technology, 1993, 115(3): 213-220. DOI:10.1115/1.2905996
[9]
邢朱苗. 大型电站凝汽器三维数值研究[D]. 保定: 华北电力大学(保定), 2009.
XING Z M. Three-dimensional numerical value research on the large power station condenser[D]. Baoding: North China Electric Power University (Baoding), 2009. (in Chinese)
[10]
PETUKHOV B S. Heat transfer and friction in turbulent pipe flow with variable physical properties[J]. Advances in Heat Transfer, 1970, 6: 503-564.
[11]
HU H G, ZHANG C. A new inundation correlation for the prediction of heat transfer in steam condensers[J]. Numerical Heat Transfer, Part A: Applications, 2008, 54(1): 34-46. DOI:10.1080/10407780802024963
[12]
GHIAASIAAN S M. Two-phase flow, boiling and condensation: In conventional and miniature systems[M]. New York: Cambridge University Press, 2008.
[13]
ZHANG C, BOKIL A. A quasi-three-dimensional approach to simulate the two-phase fluid flow and heat transfer in condensers[J]. International Journal of Heat and Mass Transfer, 1997, 40(15): 3537-3546. DOI:10.1016/S0017-9310(97)00014-8
[14]
MOUKALLED F, MANGANI L, DARWISH M. The finite volume method in computational fluid dynamics: An advanced introduction with OpenFOAM® and Matlab®[M]. Cham: Springer, 2016.