错流板翅式换热器分布参数模型构建和优化研究

李科,文键,厉彦忠,Andrea DIANI

(1.西安交通大学能源与动力工程学院,710049,西安;2.Department of Industrial Engineering, University of Padova, 35131, Padova, Italy)

板翅式换热器是一种高效、紧凑、轻巧的换热设备[1],已在航空航天、石油化工、冶金、动力工程、低温分离和液化等领域得到广泛的应用[2]。很多文献采用CFD的方法研究翅片单通道的流动换热性能[3-7]。对板翅式换热器整体进行CFD数值计算需要巨大的网格量,因此一些文献采用集总参数法[8-12]进行研究,即对于某一个通道中的流体物性计算都采用一个定性温度,该方法具有一定的局限性。在超过两股流体的板翅式换热器中,采用分布参数法更加合适[13-16],即把板翅式换热器沿流动方向分为若干小段,假定每一个离散段内工质热物性近似不变。若翅片通道中流体都是顺流或者逆流,则模型可被简化为二维。但是,在机载环境下,为了匹配各股流体的阻力,经常将冲压空气与其他工质的流动方向错流布置。文献[17-21]讨论并提出了多股流错流板翅式换热器的数学模型,将换热器沿着长度和宽度方向划分为多个单元格组,以一个单元格组为研究对象,分别对各层通道横向和纵向建立能量守恒方程。目前,错流板翅式换热器的研究很少引入轴向导热和横向导热因素,且很少涉及整体结构参数和工况参数的优化研究。

本文构建了多股流错流板翅式换热器的三维分布参数计算模型,考虑流体物性变化和隔板中轴向横向导热效应,采取高斯-赛德尔迭代法求解方程组,并将数值计算结果与热动实验台中的三股流板翅式换热器实验结果进行比较。在此基础上,研究了错流板翅式换热器整体结构参数和工况对换热量、泵功和单位泵功换热量的影响规律。

1 分布参数模型

1.1 模型简化

错流板翅式换热器的三维分布参数模型如图1所示,该模型考虑了流体物性变化及其隔板的轴向横向导热效应并假设:①流体入口处均匀分布;②稳态;③仅考虑隔板、翅片固体域及翅片通道中流体域;④未考虑翅片厚度方向和隔板厚度方向的导热;⑤翅片只考虑沿着翅高方向的导热;⑥忽略流体中导热项;⑦未考虑相变;⑧换热器与周围环境绝热。

(a)三维模型

图1(a)是错流板翅式换热器的三层翅片模型,图中将隔板和流体通道沿着x方向和y方向分割为多个小单元体。每个小单元体中均布置了储存流体温度或隔板温度的节点,在每层翅片通道的每个小单元体中,采用集总参数法,即取该小单元体中的温度和压力作为定性温度和定性压力。若换热器的体积很大,且轴向和横向的温度梯度很大,则可将换热器分割为更多的小单元体,此时仍可以保证模型的计算精度。以最下面一层翅片通道为例,如图1(b)所示,图中,实心节点是隔板中的温度节点,空心节点是流体通道中的节点,沿着流动方向隔板节点和流体节点错开半个网格长度。流道中的流体单元由两个流体温度节点(f,i,j,k)和(f,i,j+1,k)定义,隔板中的控制容积由隔板中的一个温度节点表征,如图2(a)所示。图2(b)是后续使用的一些翅片参数,图中,h是翅高,s是翅距,t是翅厚。

(a)B视角的节点布置及控制容积

图3为局部的温度节点示意图。沿着x方向(流动方向)的编号是j;沿着z方向是翅片通道层的变化,编号为i;沿着y方向是横向的导热和对流换热,编号是k。图3表明,隔板中的温度节点Tw,i,j,k与上下流道中的4个节点(Tf,i-1,j,k、Tf,i-1,j+1,k、Tf,i,j,k和Tf,i,j+1,k)、横向的相邻隔板节点(Tw,i,j,k+1和Tw,i,j,k-1)、前后的相邻隔板节点(Tw,i,j+1,k和Tw,i,j-1,k)以及上下相邻隔板的节点(Tw,i-1,j,k、Tw,i+1,j,k)有关,该隔板节点所在的控制容积的能量守恒原则可将这些温度节点联系起来。

图3 错流板翅式换热器局部温度节点示意Fig.3 Local Temperature node distribution of cross-flow plate fin heat exchanger

1.2 模型构建

流道中的流体通过隔板来换热,而流体和隔板的换热分为两部分,分别是翅片和隔板的换热(图4中第1部分,翅片和隔板之间的换热来源于翅片和流体之间的换热,这部分换热量来自二次换热面积,即翅片表面)以及流体和隔板的换热(图4中第2部分,一次换热面积)。为了求得隔板和翅片之间的导热量,需先求得翅片中的温度分布,进而根据傅里叶导热定律得到翅片根部的导热量。在实际情况中,翅片中存在沿着流体流动方向和翅高方向的导热。为了简化计算,只考虑沿着翅高方向的导热(图4中沿着z方向)。

图4 翅片一维导热和对流换热Fig.4 1-D heat conduction in fin and heat convection

翅片一维导热方程中,流体和翅片之间的换热量可看作换热源项,在温度节点(f,i,j,k)和(f,i,j+1,k)所定义的流体单元中,翅片导热方程所满足的边界条件是

z=0,T=Tw,i+1,j,k;z=h,T=Tw,i,j,k

(1)

求解给定边界条件的翅片一维导热方程,然后根据傅里叶导热定律可得翅片顶部和根部的温度梯度

(2)

式中:θi-i,j,k、θi-i+1,j,k是翅片顶部和底部的过余温度,K;me,i,j,k是求解翅片一维导热方程过程中的无量纲中间参数;θi-i,j,k、θi-i+1,j,k和me,i,j,k的计算式为

(3)

隔板所在的控制容积以及流体单元控制容积应当满足能量守恒,基于此推导得到以影响系数形式表达的求解隔板温度节点和流体温度节点的方程组。为了简化影响系数的表达式,定义

(4)

基于能量守恒的思想,推导得到计算隔板温度节点的方程组

(5)

式(5)表示温度节点Tw,i,j,k的计算依赖于上下相邻隔板、流道以及同一层隔板中(包含轴向和横向)的4个温度节点的影响。影响系数的计算方法为

(6)

式中:影响系数Aw,i,j-1,k和Aw,i,j+1,k代表了隔板中轴向相邻温度节点Tw,i,j-1,k和Tw,i,j+1,k对温度节点Tw,i,j,k的影响,相当于考虑了隔板中的轴向导热效应；影响系数Aw,i,j,k-1和Aw,i,j,k+1代表了隔板中横向相邻温度节点Tw,i,j,k-1和Tw,i,j,k+1对温度节点Tw,i,j,k的影响,相当于考虑了隔板中的横向导热效应;λw,i,j-1/2,k是隔板温度节点(w,i,j,k)和(w,i,j-1,k)的交界面上的导热系数(W/(m2·K)),是相邻两节点处导热系数的几何平均值[22];tclip是隔板的厚度,m;N是沿着x方向分割的单元体数量;w是沿y方向分割的单元体数量。

其他影响系数的计算表达式为

(7)

式(7)采用了简化的系数表达式,其中的参数bi,j,k、di,j,k、ci,j,k可自动赋予式(4)中参数λs、me、hc下标i、j、k,赋予式(4)中的参数t、h、s下标i,W、L、N、w不需要下标,λs,i,j,k是根据隔板温度节点Tw,i,j,k得到的固体导热系数。对于图4中最下面的一层翅片通道,流体沿着x方向流动,但是在横向(y方向)存在对流换热热阻和流体导热换热热阻。文献[17]考虑在错流换热器中引入了这两方面热阻的量化值,结果发现流体单元中横向能量增量相比于流动方向能量增量非常小,因此在错流板翅式换热器的计算中往往忽略了横向温度节点的影响。

流道中的流体单元(由流道中温度节点(f,i,j,k)和(f,i,j+1,k)所定义)满足能量守恒方程

(8)

式(8)经过推导,得到流体通道中节点温度的计算方法(流动方向沿着x方向)

(9)

方程组(9)中前3项是主要影响节点,表征了流动方向的能量增量,后两项表征了横向的相邻流体温度节点的影响(忽略不计)。影响系数的计算方法为

(10)

若流动沿着-x方向,也可根据能量守恒推导得到计算流体温度的方程组

(11)

式中影响系数的计算与式(10)类似

(12)

隔板中节点温度的求解采用点迭代法中常用的高斯-赛德尔迭代法[22]。为加快计算速度,迭代过程中可适当超松弛[22]。在每一个流体单元中,采用Manglik 和 Bergles[23]的经典关联式来确定流体换热系数以及摩擦阻力,第i层通道中由流体温度节点(f,i,j,k)和(f,i,j+1,k)所定义的流体单元中的压力降包含了摩擦压降以及动量变化率所产生的压降两部分,计算式[2]为

(13)

(14)

为了计算式(5)、(9)、(11)中某些影响系数,需要计算流体单元中的j因子。同样,采用文献[23]的关联式来确定锯齿翅片通道中的每个流体单元中的j因子,j因子与翅片结构参数和工况参数Re有关,Re的计算参见式(14),而根据j因子可计算该流体单元中的换热系数hc

(15)

式中λf是流体的导热系数。

错流板翅式换热器的换热量即热流体入口的焓流减去出口的焓流,或者冷流体出口的焓流减去冷流体入口的焓流。另一方面,需要考虑把单位代价所获得的换热量表示为一个无量纲量,可把换热器的流动阻力即压降转化为泵功(换热量比泵功是无量纲的)。

假设在板翅式换热器中的每一层通道中只考虑其在板束中的阻力,在某层翅片通道的管路中连接了一台泵,通过泵的做功,换热器通道出口的压力重新上升为此翅片通道的入口压力,公式[24]为

(16)

式中:Hp是泵的扬程,m;g是重力常数,9.8 N/kg;ρout和ρin是流道进、出口的密度,kg/m3。泵的有效功率与通道中的质量流量和扬程有关

Np=mgHp

(17)

式中:m是质量流量,kg/s;Np是泵功,W。

1.3 非线性方程的求解方法

迭代计算方法如图5所示,在正式进入迭代轮次之前先进行物性拟合,拟合的数据来源于nist数据库,还要进行温度场和压力场的初始化。

图5 计算流程Fig.5 Calculation process

外迭代轮次的第一个计算模块是压力场的内迭代计算(根据式(13))确定新的压力场,此时根据新的压力场再更新一次流体物性。第二个计算模块是求解方程组(5)、(9)、(11),方程组(5)的求解采用高斯-赛德尔迭代法,方程组(9)、(11)是步进计算。每一个外迭代轮次结束之后,进行温度残差的判断,残差满足收敛条件即可输出结果。反之,则进入下一个外迭代轮次,根据更新后的温度场和压力场确定各流体单元中的物性。

1.4 实验验证

参阅文献[25]所给的三股流板翅式换热器实验装置,芯体长0.4 m(沿着x方向)、宽0.13 m(沿着y方向),工况参数见表 1,采用了ABCABCABC……ABCA的排布方式),翅片结构参数可参见文献[25]。在实验中,改变表1中A流体的入口流量,在不同的流量下测定B和C流体的出口温度,将程序的计算结果和实验结果[25]进行比较,如图6所示。

表1 错流板翅式换热器工况参数

图6表明,热流体B和C的出口温度随着冷流体A流量的增大而增大,理论计算表明热流体B和C的出口温度的差异小于实验结果。这是由于B和C流体在流动过程中也会互相换热,温度趋于一致。理论假设流道中的流体均匀分布,但实际中B和C流体通过封头和导流片进入换热器后的分布仍然不均匀,导致B和C流体未经过充分换热。因此,图6中计算值和实验值具有0.2～3.0 K的温差是合理的。

图6 错流板翅式换热器实验值[25]和模拟值的比较Fig.6 Numerical results versus experiment results in cross-flow plate fin heat exchanger

2 计算结果分析

采用1.4小节实验验证所采用的错流板翅式换热器进行讨论。图7是A流体不同流量和B流体不同入口温度下板翅式换热器的换热量的变化,图中TBin为B流体入口温度。可以看出,换热量随着A流体流量的增大而增大,但是增大的速度趋缓。A流体的流量增大,A侧的对流换热系数增大,则A流体作为冷流体能够吸收更多的热量,因此换热量增大。当TBin取323～403 K时,换热量随着A流体流量增大而增大的趋势类似。当TBin为303 K时,换热量随着A流体流量增大而增大的速度明显小于当TBin取其他温度时的。例如,当A流体的流量从800 kg/h增加至2 400 kg/h时,TBin取303 K,换热量增大了13.3%,而TBin取323 K,换热量增大了31.9%。变化趋势不同的原因是:若TBin为303 K,与A流体的入口温度(301.6 K)非常接近,A流体与B流体几乎不发生热交换,主要是A流体和B流体从C流体中获取热量,尽管A流体流量从800 kg/h增加至2 400 kg/h,增加量很大,但是C流体的流量很小,成为了限制换热量增加的瓶颈。

图7 换热量随A流体流量和B流体入口温度的变化Fig.7 Variation of heat transfer amount with mass flowrate of fluid A and inlet temperature of fluid B

在换热器实际运行过程中,还需考虑其泵功消耗。图8是换热器在不同A流体流量和不同TBin下泵功消耗的变化。可以看出,不同TBin情况下,泵功消耗随着A流体流量增大而增大的趋势相同。A流体流量增大,表明A流体的流速增大,导致流动阻力增大,必然导致泵功消耗的增大。在A流体流量相同的情况下,增加TBin,则B流体的密度减小,流速增加,使得B流体通道中的流动阻力增大,即换热器的泵功消耗增大。

单位泵功换热量也是一个值得关注的评价指标。单位泵功换热量随A流体流量和B流体入口温度的变化如图9所示。可以看出,当A流体的流量不变时,TBin从323 K增大至403 K,单位泵功换热量一直增大。这从侧面反映出换热量随着TBin增大而增大的速度要大于泵功消耗随着TBin增大而增大的速度。若保持TBin不变,单位泵功的换热量先随着A流体流量的增大而增大。当A流体的流量增大至1 300 kg/h时,单位泵功换热量达到最大值,随后单位泵功换热量随着A流体的流量的增大而减小。综合分析可见,选择合理的运行工况能提高板翅式换热器的工作效率。

前面讨论的板翅式换热器的芯体长度是0.4 m,宽是0.13 m,长(沿x方向)和宽(沿y方向)的乘积是0.052 m2,这个面积是板翅式换热器的芯体在xy平面上的投影面积。保持xy面投影面积不变,可研究当y边长度在0.1～0.5 m范围内变化时换热器换热量、泵功消耗和单位泵功换热量的变化,然后讨论在不同的xy面投影面积下,换热器单位泵功消耗量随着换热器y边长度的变化。

图10是换热器在不同xy面投影面积下,换热量随着换热器y边长度的变化。可以看出,当换热器xy面投影面积不变时,换热器的换热量基本不发生变化。例如,当换热器xy面投影面积是0.06 m2时,换热器y边长度从0.1 m增大至0.6 m,换热量的变化范围大约是此时换热量的3.21%,这种变化更可能是计算误差所引起的波动。当控制换热器y边长度为0.25 m时,换热器xy面投影面积从0.04 m2增大至0.07 m2,换热量增大了7.3%,即对于此错流换热器,换热器面积增大75%(当换热器沿着z方向的高度不变时,xy面投影面积实际上表征了总的换热面积)仅仅使得换热量提升了7.3%,提升效果很不显著。1.4小节实验验证所采用的错流换热器的xy面投影面积是0.052 m2,但是图10中的数据表明:当y边长度是0.25 m时,xy面投影面积为0.04 m2的换热器的换热量相对于xy面投影面积为0.052 m2的换热器的换热量仅下降了3.81%。因此,在保持换热器沿着z方向高度不变且只关注换热量的情况下,把xy面投影面积设计为0.04 m2也未尝不可。

图10 不同xy面投影面积下换热量随着y边长度的变化Fig.10 Variation of heat transfer amount with y direction length under different xy projected area

文中讨论的y边长度方向是A流体的流动方向,A流体的流量要比B和C两股流体的流量之和都要大。当保持换热面积不变时,增大A流动方向的长度意味着A的长度增大且A的横向的宽度减小,不仅导致A流体的流速增大,且A流体要流过更长的距离,结果导致整体的泵功消耗急剧增大,如图11所示。这表明在设计错流换热器时,对大流量流体流动方向的长度不能设计得过长,否则为了减小该股流体的泵功消耗,势必要增大其横向宽度,导致换热面积增大。例如,在图11中控制y边长度为0.35 m,泵功的消耗随着xy面投影面积的减小而增大。图11还表明,当A流体流动方向(y边长度)的长度取0.15～0.20 m时,泵功的消耗最低。

图11 不同xy面投影面积下泵功消耗随y边长度的变化 Fig.11 Variation of pump power with y direction length under condition of different xy projected area

图12是错流板翅式换热器单位泵功换热量的变化,即换热器换热量与泵功的比值的变化。由于图10中的换热量随y边的长度基本不发生变化,图11中的泵功随着y边长度的变化剧烈,导致单位泵功换热量的变化较为突兀。随着板翅式换热器y边长度的增大,单位泵功换热量先增大后减小。当y边的长度取0.15～0.20 m时,单位泵功换热量达到最大值。实际上,在整个y边长度变化范围内可以认为换热量不发生变化,主要是y边长度的增大导致了泵功消耗急剧增大,使得单位泵功换热量急剧减小。可见,当控制错流板翅式换热器的xy面投影面积不变时,适当选择大流量流体沿流动方向的换热器长度,可以提高单位泵功换热量。

图12 不同xy面投影面积下单位泵功换热量随y边长度的变化Fig.12 Variation of heat transfer amount per unit pump power with y direction length under condition of different xy projected area

3 结 论

本文构建了多股流错流板翅式换热器的分布参数模型,考虑隔板中的轴向和横向导热效应以及流体通道中的物性变化,采用MATLAB编写了程序,讨论了其换热量、泵功和单位泵功换热量的变化,获得了使泵功最小及单位泵功换热量最大的整体结构参数和工况参数,本文结论如下。

(1)所研究的错流板翅式换热器的换热量和泵功消耗都随A流体的流量和B流体的入口温度增大而增大,当A流体增大至一定程度,换热量的增大趋缓。

(2)在所讨论的错流板翅式换热器模型中,单位泵功换热量存在极值,先随着A流体流量的增大而增大,当A流体的流量增大至1 300 kg/h时,又随着A流体流量的增大而减小。

(3)对错流板翅式换热器,维持其xy面投影面积不变,那么由于换热面积不发生变化,整体的换热量基本不变,但是泵功消耗会随工况发生变化。计算结果表明,采用4种不同的xy面投影面积,当换热器y边长度为0.15～0.20 m时,泵功消耗最低,单位泵功换热量最大。