中冷器作为一种典型的换热器, 其作用是冷却来自增压器的高温高压空气，采用增压中冷的方式可以大幅度降低污染物的排放并提升柴油机的动力性能^[1-2].当前，关于中冷器的研究主要集中在如何能在确保其对增压空气的降温效果的同时降低其压降.为了节约成本，缩短研发时间，目前多数研究者采用fluent等多维仿真软件对换热器的换热量与压力损失进行设计计算，如Chen等^[3-5]通过对不同种类翅片进行试验拟合的方式得出数种不同结构翅片的传热因子与摩擦因数关于雷诺数的关联式，可以较好地评价翅片的换热效率, 为换热器设计提供了便利; 张钦国^[6]对某型号中冷器的CFD流场仿真时发现，在弯管入口区域存在回流，且回流随来流增加而增大，这是入口压力损失的主要来源之一.多数研究者^[7-9]将换热器CFD仿真值与试验值的偏差归结为随机误差或由于模型过度简化所带来的偏差，但笔者认为基于压力求解的方式忽略了密度的变化，即不考虑流体的可压缩性, 这也是不容忽视的偏差来源之一.

本文以某轮式装载机水冷中冷器增压空气侧为研究对象, 建立了基于翅片壁面摩擦，换热和流动截面变化的换热器可压缩流模型; 将UDF文件导入fluent进行三维流场仿真, 并与同条件下基于压力求解的仿真结果进行对比, 验证模型的可靠性; 将可压缩模型与不可压缩模型的计算结果与某换热器厂家的中冷器台架试验数据进行对比，研究可压缩流动对换热器计算的影响.

1 换热器可压缩流模型

影响管内气体流动的因素有很多，如壁面摩擦、重力引起的质量力、质量变换、热交换等，但在可压缩流动中并不是每种影响因素都起着关键性作用, 而是根据工况不同影响因素有主有次^[10].流体介质为空气的换热器，在空气流经换热器时通常会产生较大的换热量和压力损失，因此根据换热器的特性,本文以空气侧摩擦阻力、空气侧换热量、换热器流道截面变化为研究对象,建立换热器可压缩流模型.

1.1 存在热交换的流动模型

流体运动的基本守恒定律参见文献[11]，本节着重分析只存在换热时的可压缩流对换热器计算的影响.建立存在换热的可压缩流模型，在一维定常流假设下对中冷器模型进行简化，认为其内部流动是无机械功变化的可压缩无摩擦换热管流，并引入式(1)~(4)^[12]：

$ \tau \left( \lambda \right) = 1 - \frac{{k - 1}}{{k + 1}}{\lambda ^2}, $

(1)

$ \pi \left( \lambda \right) = {\left( {1 - \frac{{k - 1}}{{k + 1}}{\lambda ^2}} \right)^{\frac{k}{{k - 1}}}}, $

(2)

$ q\left( \lambda \right) = {\left( {\frac{{k + 1}}{2}} \right)^{\frac{1}{{k - 1}}}}\lambda {\left( {1 - \frac{{k - 1}}{{k + 1}}{\lambda ^2}} \right)^{\frac{k}{{k - 1}}}} = \frac{{{A_{{\rm{cr}}}}}}{A}, $

(3)

$ f\left( \lambda \right) = {\left( {\frac{2}{{k + 1}}} \right)^{\frac{1}{{k - 1}}}}q\left( \lambda \right)z\left( \lambda \right). $

(4)

式中：冲量函数z(λ)=λ+1/λ，k为绝热指数，A_cr为临界截面面积，λ为纲量一的速度系数

由能量方程与动量方程联立得到式(5)的滞止温度比：

$ \frac{{T_b^ * }}{{T_a^ * }} = \left[ {\frac{{z\left( {{\lambda _a}} \right)}}{{z\left( {{\lambda _b}} \right)}}} \right], $

(5)

将由连续方程求解得到的密度比代入得到温度比

$ \frac{{{T_a}}}{{{T_b}}} = \frac{{T_b^ * \tau \left( {{\lambda _b}} \right)}}{{T_a^ * \tau \left( {{\lambda _a}} \right)}} = {\left[ {\frac{{z\left( {{\lambda _a}} \right)}}{{z\left( {{\lambda _b}} \right)}}} \right]^2}\frac{{\tau \left( {{\lambda _b}} \right)}}{{\tau \left( {{\lambda _a}} \right)}}, $

(6)

由(4)(7)联立可以得到截面a、b的滞止压力比与压强比:

$ \frac{{p_b^ * }}{{p_a^ * }} = \frac{{f\left( {{\lambda _a}} \right)}}{{f\left( {{\lambda _b}} \right)}} = \frac{{{p_b}\pi \left( {{\lambda _b}} \right)}}{{{p_a}\pi \left( {{\lambda _a}} \right)}}. $

(7)

式中：p^*为滞止压力，T^*为滞止温度.

将上述模型与基于压力(不可压缩模型)的求解值进行对比.在基于压力求解时, 换热器的压力损失Δp与换热量Q分别为

$ \Delta p = \frac{1}{2}\rho fU_{\min }^2, $

(8)

$ Q = {c_p}\dot m\Delta t. $

(9)

式中：f为摩擦因数，U_min为来流速度与流通率之比，c_p为比定压热容，${\dot{m}}$ 为质量流量，Δt为温差.

不可压缩模型的摩擦因数f为Kays & London的矩形错齿翅片试验关联式^[13].

在如下工况对比上述模型：某轮式装载机采用水冷中冷器冷却增压空气，视其冷侧为理想冷侧，其热侧增压空气温度为111.5 ℃，表压200 kPa，亚音速，绝热指数为1.4，芯体流通率为0.438，热侧换热面积为2.06 m².采用可压缩模型在MATLAB中求解不同进气量下的中冷器增压空气出口压力波动与换热量，计算结果如图 1所示.可见压力变化值很小且呈现不规则波动，基本可以认为中冷器在可压缩流计算中由热交换引起的压力变化很小，因此下文中不考虑由其引起的压力变化; 图 1中的换热量为可压缩模型的换热量与由式(9)计算的换热量差值，其随流量增大呈现加速上升趋势，即流量越大可压缩模型与不可压模型产生的换热量差值越大.可见，图 1中的计算结果符合Rayleigh曲线亚音速冷却阶段压力、密度及温度的变化规律.

1.2 摩擦和截面变化对压力损失的影响

可压缩流压力损失主要来自两部分：一部分为截面变化引起的损失，另一部分为壁面摩擦引起的压力损失.首先讨论由摩擦引起的压力损失.假设换热器的热侧为一维定常有摩擦绝热管流，流动过程无机械能变化.以量纲一速度系数表示平均摩擦因数f(如式(10))，结合Fanno曲线及以摩擦为自变量的压力微分式(11)，可知无论来流是超音速或亚音速, 管出口静压和全压都会呈下降趋势.因此，在进口条件已知时存在以下积分解：

$ \left( {\frac{1}{{\lambda _a^2}} - \frac{1}{{\lambda _b^2}}} \right) - \ln \frac{{\lambda _b^2}}{{\lambda _a^2}} = \frac{{8k}}{{k + 1}}\bar f\left( {\frac{L}{D}} \right), $

(10)

$ \frac{{{\rm{d}}p}}{p} = - \frac{{kM{a^2}\left[ {1 + \left( {k - 1} \right)M{a^2}} \right]}}{{2\left( {1 - M{a^2}} \right)}}4f\frac{{{\rm{d}}x}}{D}. $

(11)

使用临界截面cr表示平均摩擦系数和总压力比，可压缩流由摩擦引起的压力损失可以按式(12)和(13)计算:

$ \frac{p}{{{p_{{\rm{cr}}}}}} = \frac{1}{{Ma}}{\left[ {\frac{{k + 1}}{{2\left( {1 + 0.5\left( {k - 1} \right)M{a^2}} \right)}}} \right]^{1/2}}, $

(12)

$ \frac{{{p_{\rm{b}}}}}{{{p_{\rm{a}}}}} = {\left( {\frac{p}{{{p_{{\rm{cr}}}}}}} \right)_b}/{\left( {\frac{p}{{{p_{{\rm{cr}}}}}}} \right)_a}. $

(13)

在考虑截面变化带来的出口压力变化时, 可将换热器热侧简化为类似拉伐尔管的模型，以上文中冷器为例，由于已知空气压缩器的出口压力与流量，将中冷器流通率代入流量函数q(λ), 结合式(3)便可直接求解出口与入口的压力比, 进而得出截面变化对压力的影响.代入1.1节的工况条件, 结合该水冷中冷器的几何尺寸, 计算出进气压力为200 kPa时压损约为4.25%.将修正的两部分压力损失叠加, 在MATLAB中计算压力损失并与按式(8)的计算值进行对比, 结果如图 2所示，图 2中Kays & London曲线的压力损失按文献[13]中3/4-11.1型锯齿翅片摩擦因数计算.可知在来流速度较低时, 流体的可压缩性并不会对热侧压力损失产生明显影响; 但在来流速度超过20 m/s时, 两者偏差值已超过了5%.由此可见，流体力学中认为来流速度低于0.3 Ma不必考虑流体的可压缩性并不能很好地适用于流体介质为空气的换热器计算.

1.3 换热器的可压缩流模型

考虑多重影响因素共同作用下的可压缩流模型，对于完全气体的流动，忽略重力的影响可以将基本守恒定律改写成如下形式：

能量方程可以改写为

$ \delta Q - \delta W - {\rm{d}}{H_{\rm{a}}} = {\rm{d}}H = {c_p}{\rm{d}}{T^ * }. $

(14)

式中：H=h+v²/2+gz，dH_a=(H－H_a)dm/${\dot{m}}$，z为重力g方向的位移.将式(14)除以c_pT得

$ \frac{{\delta Q - \delta W - {\rm{d}}{H_a}}}{{{c_p}T}} = \left( {1 + \frac{{k - 1}}{2}M{a^2}} \right)\frac{{{\rm{d}}{T^ * }}}{T}. $

(15)

动量方程可以改写成

$ \frac{{{\rm{d}}p}}{p} + kM{a^2}\frac{{{\rm{d}}Ma}}{{Ma}} + \frac{{kM{a^2}}}{2}\left( {\frac{{{\rm{d}}T}}{T} + \frac{{4fdx}}{D}} \right) + kM{a^2}\left( {1 - y} \right)\frac{{{\rm{d}}\dot m}}{{\dot m}}. $

(16)

状态方程为

$ \frac{{{\rm{d}}p}}{p}{\rm{ = }}\frac{{{\rm{d}}\rho }}{\rho } + \frac{{{\rm{d}}T}}{T}. $

(17)

连续方程为

$ \frac{{{\rm{d}}\dot m}}{{\dot m}} = \frac{{{\rm{d}}\rho }}{\rho } + \frac{{{\rm{d}}V}}{V} + \frac{{{\rm{d}}A}}{A}. $

(18)

联立式(15)~(18)可解得马赫数Ma与温度T、截面变化dA、水力直径D、摩擦因数f等流动参数之间的关系式：

$ \begin{array}{*{20}{c}} {\frac{{{\rm{d}}Ma}}{{Ma}} = \frac{{\left( {1 + \frac{{k - 1}}{2}M{a^2}} \right)}}{{1 - M{a^2}}}\left[ {\frac{{kM{a^2}}}{2}\left( {\frac{{4f{\rm{d}}x}}{D}} \right) - } \right.}\\ {\left. {\frac{{{\rm{d}}A}}{A} + \frac{{1 + kM{a^2}}}{2}\frac{{{\rm{d}}{T^ * }}}{{{T^ * }}}} \right].} \end{array} $

在多数工程问题中，换热器的几何参数与进口参数为已知条件，设定计算步长为Δx，本文用Runge-Kutta迭代法求解延长度方向任意截面上的马赫数，再带入式(6)、(7)中可求解多重影响因素共同作用下的各截面的其它流动参数.

2 仿真与试验分析 2.1 仿真模型及边界条件

采用fluent软件对某中冷器热侧进行可压缩多维流动分析.将前文讨论的多重影响因素下的压力与温度求模型编写的UDF文件导入fluent中进行三维流场仿真.

本文使用的某装载机水冷中冷器热侧由10层散热带组成，冷侧由9层热板组成.中冷器的实体模型如图 3所示.

其芯体尺寸为120 mm×300 mm×110 mm，冷热流体流动形式见图 4.增压空气沿水平方向流动，冷却液为二次U形流动，中冷器热侧矩形错齿翅片节距为2.4 mm，翅片高6.5 mm，冷侧矩形错齿翅片节距为2.5 mm，翅片高3 mm.使用Gmbit对其进行网格划分，由于使用多孔介质模型模拟芯体部分对网格要求不高，故将模型网格精度控制在0.8以下即可，按文献[14]进行网格无关性验证，模型最终网格数量为189万.分别设置求解器为基于压力和基于密度求解，基于压力求解的流动起始点增压空气物性参数如表 1所示^[15].

如图 5、6所示，在使用多孔介质模型时，基于密度求解与基于压力求解在温度场上并不存在明显的差异，两者的出口温度几乎一致，这是因为使用多孔介质模型时换热器芯体失去了内部结构的概念，只按设定的孔隙率、黏性阻力因子和惯阻性阻力因子进行模拟计算，这就造成了多孔介质模型对可压缩流动的“不敏感”.笔者认为这是一些使用多孔介质模型来模拟换热器的文献往往出现仿真值和试验结果存在较大偏差的原因之一.

因此, 若要较好地反映可压缩流的影响, 必须考虑换热器内部翅片结构; 但若在保证网格质量的前提下, 对包含翅片的10层热侧散热带划分网格则网格数量将高达上千万，显然目前的计算机硬件水平难以实现对这个量级的模型进行求解^[16].本文以中冷器热侧的一层散热带建立如图 7所示的简化模型，其网格质量为0.35，网格数量为109万.为了防止中冷器进出口突扩造成的扰动, 该模型省略了中冷器的进、出口突缩与突扩段，将冷侧隔板视为理想冷侧.在不考虑每层散热带之间干涉的情况下, 热侧每一层散热带的压损将与中冷器进出口压损一致.导入前文UDF文件进行可压缩流仿真计算，作为对照组的不可压缩流依然使用前文的进口参数.

2.2 仿真结果分析

图 8为散热带中心截面温度云图，图 9和图 10分别为压力和换热量曲线.

由图 8可知, 换热量随流量增加而增加，但增加速度呈现逐渐减缓趋势.这是因为热侧传热效率随着温度升高而降低，当进气量持续增加时, 换热量增量将逐渐趋于平缓.与前文的一维解析解不同, fluent中两模型的换热量差值远大于前文的一维解，由图 1可知, 当进气量达到0.5 kg/s时, 可压与不可压换热量差值应在0.6 kW左右，而如图 9所示的fluent仿真结果表明，换热差在同流量下达到了2 kW左右.笔者认为其原因主要有3点：首先，一维流动旨在研究流动参数沿轴线上的变化，而fluent计算则涵盖了流体沿各方面的变化，这使得三维仿真解的精度与完整程度要高于一维计算; 第二，对芯体的拉伐尔管简化, 忽略了翅片间流场的相互影响，这对换热量计算有一定的影响; 第三，压力损失仿真结果如图 10所示，两条曲线压损随着流量的增加呈加速上涨状态，但考虑压缩性时压力损失的增长速度低于传统的不可压缩曲线，其主要原因是密度的增加造成动压下降，因而减少了摩擦损失，若不考虑压缩性，则压力损失与速度的平方成正比.

2.3 试验对比分析

委托某换热器厂商进行换热器台架试验.试验标准参照JB/T 8577-2015《内燃机水散热器技术条件》，在中冷器热侧进、出口安装压力与温度传感器，冷侧进、出口安装温度传感器，用DEWE-43和krypton-TH8数据采集仪采集数据，使用DEWESoft-X2软件进行数据后处理.试验条件如下：热侧进气状态为200 kPa，110 ℃，冷侧乙二醇冷却液流量为100 L/min，进温为40 ℃.

可压模型与不可压缩模型的换热量与压力损失CFD仿真与试验值对比见图 11、12所示，压力和换热量对比试验值的偏差见表 2，其计算公式为

$ 偏差 = \frac{{仿真值-试验值}}{{试验值}} \times 100\% . $

当空气流量超过0.2 kg/s时, 使用不可压缩模型CFD仿真值与试验的值偏差开始持续增大，而可压缩曲线依然能与试验值保持较高的吻合度，验证了前文中提出的中冷器可压缩模型的准确性.通过换算得到:在流量为0.2 kg/s时, 速度约为0.15 Ma; 在0.3 kg/s时速度为0.19 Ma; 在速度低于0.1 Ma时, 两条仿真曲线与试验值的偏差几乎可以忽略; 当速度达到0.15 Ma时, 不可压缩模型与试验值的偏差为5%;而当速度达到0.2 Ma时, 不可压缩模型换热量和压力损失与试验值偏差值都超过了10%;随着速度继续增加至超过0.3 Ma时, 压力损失偏差值甚至超过了30%;而实际工程问题的最高空气流速多在0.2~0.3 Ma，意味着使用不可压缩模型计算时会产生一定的偏差.笔者认为本文的CFD仿真与试验值差异主要来自以下3方面：首先，不可压缩模型使用的矩形错齿翅片壁面摩擦因子为Kays & London的3/4-11.1型锯齿翅片试验关联式，其几何参数较为接近本文研究的中冷器热侧翅片, 但并不完全相同.在低雷诺数区域能很好地吻合试验值，但在高雷诺数区域则与试验值产生了较大的偏差，这就造成了图 11在高流量时的不可压CFD值与试验值差距不断增大; 其次, CFD仿真忽略了中冷器进出口突缩突扩及相邻两层散热带之间的干涉影响; 最后, CFD仿真时对冷侧的理想化假设造成了换热量与试验值的差异，试验冷却液流量为100 L/min, 在高进气量下势必会造成换热效率下降，导致试验换热量低于仿真换热量.

3 结论

1) 基于壁面摩擦、换热、流道截面变化建立了换热器可压缩流模型，与传统基于不可压求解方式相比，可压缩流模型与试验值匹配精度更高，其中换热量与试验值最大偏差为6.5%，压力损失最大偏差为7.5%，为换热器优化设计提供了一定指导.

2) 在进气流量较小时, 基于不可压求解能较好地吻合试验值，但随流量的增加, 其偏差值逐渐增大; 而基于可压缩流求解, 则能在较高流速下依然保持较高的求解精度.

3) 在流体介质为空气时, 来流速度在0.15 Ma以下, 不可压模型可以代替可压模型来简化计算，提高效率; 而当来流＞0.15 Ma，就有必要考虑空气的可压缩性, 以提高计算精度.