单侧加热方形通道内超临界水传热研究

引 言

再生冷却一般利用燃料的热沉冷却燃烧室壁面，然而随着飞行速度的增加，高超声速飞行器燃烧室壁面的散热面临极大挑战。常规的解决方法是通过热管理提高燃料热沉的利用效率[1]，然而仅依靠提高燃料热沉的利用效率无法满足冷却要求，需要增加冷却剂以辅助燃料冷却燃烧室。常用冷却剂有正癸烷、航空煤油RP-3、水和其他高密度烃。与其他冷却剂相比，水具有成本低、来源广和换热能力强的特点，作为优质冷却剂广泛应用于各种换热器。由于燃烧室壁温高，常压下的水容易出现相变，而超临界水的密度近似液体的密度量级，同时又具有与气体相近的扩散能力，具有非常好的流动性和传输性能，超临界水十分适合作为飞行器燃烧室壁面冷却剂。

超临界水在核燃料堆冷却换热研究较为深入，Shang等[2]数值模拟研究了直径对超临界水在水平圆管中流动的影响，发现直径影响传热是由于边界层浮力的变化，小管径更有利于传热。范辰浩等[3]实验研究了小管径内超临界水的传热恶化特性，发现在高热通量时传热恶化的产生是由于流动层流化的发展，而管内径向的工质温度和物性的巨大差异促进了流动层流化的发展。Wang等[4]基于大涡模拟研究了圆形管道中超临界水在无重力、竖直向上和竖直向下流动下的传热，发现与无重力流动相比，向下流动的传热略有改善，向上流动的传热出现恶化，提出由于浮力的影响而产生这些差异。

Bai等[5]数值模拟了竖直向上圆形管道中超临界水非均匀加热对流动换热的影响，发现与均匀加热方式相比，非均匀加热方式表现出显著的差异，截面温度分布沿周向呈现出较大的不均匀性，传热强化仅出现在局部区域。Gao等[6]对超临界水在非均匀加热的水平圆管中的流动换热进行了数值模拟，分析了超临界水的非均匀传热特性及机理，并讨论了侧向二次流和浮升力效应对传热作用开始的判据，发现二次流与热物性之间的相互影响对传热具有重要影响。许多文献对超临界水传热强化和恶化产生的机理进行了分析[7-11]，发现流体物性变化对超临界水传热具有重要影响。文献中对超临界水换热关联式的适用性也进行了研究[12-17]，发现超临界水的关联式需要根据不同工况的换热特点进行修正。

现有研究主要关注超临界水在圆形通道内的流动传热，非圆结构如类矩形、类三角形和方环形等结构也存在相关研究[18-21]，而高超声速飞行器燃烧室壁面的冷却通道具有方形单侧加热特点，目前缺乏超临界水在方形单侧受热通道内流动传热的相关研究，亟需获得其流动传热规律。

本文对超临界水在方形冷却通道的流动换热进行数值模拟研究，通过分析流场和超临界水热物性变化，阐明传热恶化产生和恢复的机理，并验证超临界流体常用传热关联式在单侧加热方形通道内的适用性，最后，提出非对称凹槽和双通道强化结构，防止通道传热恶化，提升通道综合换热性能并进行场协同原理分析。

1 物理模型及数值模拟方法

1.1 计算模型

航空发动机冷却通道紧贴高温壁面，冷却介质在多个微小通道内流动，考虑模型的周期性，研究的模型可以简化为单侧加热的方通道，如图1所示。方通道当量直径为D，壁厚a=2 mm，方形通道加热长度为Lt，为得到较均匀的速度分布并防止回流，模拟过程中增加长度Ld=Lu=0.1 m的入口段和出口段。

图1

图1   方形通道示意图

Fig.1   Schematic of square channel

为提升通道换热性能，提出了多种强化通道。图2（a）~（c）分别表示对称凹槽、倒角凹槽和非对称倒角凹槽的结构。其中加热长度Lt=0.5 m，凹槽半径R=2 mm，槽深d=1 mm，凹槽间距为25 mm，凹槽数目为16。考虑入口处换热较强，无须强化，第一个凹槽布置在距离入口50 mm处，其余凹槽均匀分布于内壁面靠近加热面处。凹槽加倒角结构如图2（b）、（c）所示，其中圆角半径r=1 mm，r1=2 mm，r2=0.6 mm。非对称倒角结构由曲线与6个点（P0，P1，P2，P3，P4，P5）确定。此外，也研究了非对称凹槽最低点偏下游以及偏上游对通道流动换热的影响。

图2

图2   凹槽结构示意图

Fig.2   Schematic of grooves

双通道结构的布置是以1 mm厚度的壁面平均分割流体域，考虑两种布置形式，即上下布置（下层靠近加热端）与并排布置。二者均为对称结构，且结构参数一致。

1.2 边界条件与求解方法

采用质量流量入口、压力出口与均匀热流边界条件，计算工况如表1所示。基于商业软件Fluent求解，速度与压力的耦合采用SIMPLE算法求解，控制方程中对流项的离散采用二阶迎风格式，扩散项采用QUICK格式。当连续性方程、动量方程和湍流方程残差值小于1.0×10-5，能量方程的残差值小于1.0×10-8，且不随计算步数变化时认为结果收敛。超临界水热物性通过调用NIST数据库获得。

表1   计算工况参数

Table 1  Conditions of numerical simulation

1.3 模型验证

在靠近固体和流体区域的内壁附近采用精细网格，确保无量纲距离y+小于1，以满足SST k-ω模型的要求。经网格的独立性考核，最终采用627万的网格数量进行计算。

由于目前缺乏对方形通道中超临界水传热的实验研究，首先对超临界水在圆形管道中竖直向上流动的数值模拟结果与Yamagata等[22]的实验结果进行对比，实验段由AISI 316型不锈钢管光管制成，内径为7.5 mm，壁厚为2 mm。实验工况为：压力P=24.5 MPa，质量流速G=1260 kg·m-2·s-1，热通量q=233 kW·m-2。如图3所示，计算结果与文献中的数据吻合很好，最大误差为0.59%，表明本计算模型能够较好地处理超临界水的流动传热现象。

图3

图3   实验结果与数值模拟结果的比较

Fig.3   Comparison between available experimental and present numerical results

2 结果分析与讨论

2.1 传热恶化机理分析

本节对超临界水在竖直向上流动时出现传热恶化的工况（case 1）进行了分析。由于研究对象为超临界水的竖直向上流动，而此时流动方向与重力方向相反，二次流仅仅依靠流体密度差产生，影响十分微小，因此不考虑二次流动的影响。

Koshizuka等[23]对传热恶化时的表面传热系数做出如下定义:

h/hDB<0.3$h/h_{DB}<0.3$(1)

其中hDB为Dittus-Boelter关联式[24]计算的表面传热系数。

图4表示了数值模拟与Dittus-Boelter关联式计算的内壁面局部表面传热系数的比值随流体平均焓的变化，可以发现传热恶化的起始点A、传热恶化点B和传热恶化恢复点C，这是后续分析的三个基准点。

图4

图4   竖直向上与无重力流动内壁面h/hDB沿程变化

Fig.4   h/hDB on the inner wall with and without gravity

需要注意的是，当流体的焓值过高时，会再次出现h/hDB < 0.3的现象，这是由于流体温度靠近拟临界温度，比定压热容增加剧烈，Prandtl数过大，导致hDB偏高，因此前面ABC三点才是传热恶化产生的区域。

图5、图6展示了传热恶化区域A、B和C点的流场及物性分布，分别从边界层厚度与近壁区湍动能两个方面来解释传热恶化产生和恢复的原因。为防止局部高温出现，主要关注单侧加热下的受热面，故图5、图6主要表示了加热面单侧的一半流场与物性分布。其中，A到B的过程被定义为传热恶化加剧的过程，由图5（a）、（b）可以看出，从A到B点黏性底层的急剧升温导致流体动力黏度增加，动力黏度增加导致边界层中黏性底层增厚，如图5（c）所示。由于黏性底层中传热方式以导热为主，而在传热恶化区，与缓冲层与对数层相比，黏性底层内热导率整体处于较低的水平[25]，综合影响下增加了黏性底层整体热阻，降低了边界层流体的流动换热能力。而对于B到C过程，黏性底层内流体动力黏度降低，边界层厚度恢复到与A点相似的水平，从而使C点换热能力提升，传热恶化恢复到正常传热状态。

图5

图5   边界层厚度的影响

Fig.5   Effect of boundary layer thickness

图6

图6   近壁区湍动能的影响

Fig.6   Effect of near-wall turbulent kinetic energy

图6（a）可以看到，相比于A、C两点，B点处黏性底层与缓冲层内密度均处于较低水平，缓冲层内尤为明显，因此从A到B点的过程中，低密度的流体在浮力作用下加速上升，缓冲层内出现明显的热加速现象，该现象使得B点处缓冲层与对数层内的速度分布比较平缓，从而使得该处湍动能明显降低，导致传热恶化发生，如图6（b）、（c）所示。而从B到C点过程中，对数层流体密度进一步减小，流体速度增加，形成M形速度分布，增强了湍流的强度，导致湍动能增加，从而使传热恶化恢复到正常传热水平。

综上所述，导致传热恶化的原因为以下两点：（1）黏性底层内黏度增大使得边界层局部增厚； （2）黏性底层与缓冲层内密度减小导致湍动能降低。上述耦合影响导致了传热恶化的出现。

2.2 常用关联式的验证

文献中存在较多超临界流体换热的关联式，主要是在传统Dittus-Boelter关联式的形式上引入壁面物性参数修正，以提高关联式精度。选择了6个常用于超临界流体换热的关联式以验证其对方通道单侧加热超临界水的适用性。式（2）~式（7）分别来自于Gupta等[26]、Mokry等[27]、Jackson[28]、Bishop等[29]、Kim等[30]和Fan等[31]，关联式中计算的物性是根据模拟得到的流体平均温度与相应固体壁面温度，由NIST数据库查出对应的物性参数。

Nub=0.01Re0.89bPrb¯¯¯¯¯¯−0.14(ρwρb)0.93(λwλb)0.22(μwμb)−1.13$N{u}_b=0.01R{e}_b^{0.89}{\overline{P{r}_b}}^{-0.14}{\left(\frac{{\rho }_w}{{\rho }_b}\right)}^{0.93}{\left(\frac{{\lambda }_w}{{\lambda }_b}\right)}^{0.22}{\left(\frac{{\mu }_w}{{\mu }_b}\right)}^{-1.13}$(2)Nub=0.0061Re0.904bPrb¯¯¯¯¯¯0.684(ρwρb)0.564$N{u}_b=0.0061R{e}_b^{0.904}{\overline{P{r}_b}}^{0.684}{\left(\frac{{\rho }_w}{{\rho }_b}\right)}^{0.564}$(3)Nub=0.0183Re0.82bPrb0.5(ρwρb)0.3(c¯pcp,b)nn=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪0.4,TwTpc≤1或TbTpc≥1.20.4+0.2(TwTpc−1),Tb≤Tpc≤Tw0.4+0.2(TwTpc−1),1≤TbTpc≤1.2$\begin{array}{c}N{u}_b=0.0183R{e}_b^{0.82}P{r_b}^{0.5}{\left(\frac{{\rho }_w}{{\rho }_b}\right)}^{0.3}{\left(\frac{{\overline{c}}_p}{c_{p,b}}\right)}^n\\ n=\left(\begin{array}{c}0.4,\frac{T_w}{T_{pc}}\le 1或\frac{T_b}{T_{pc}}\ge 1.2\\ 0.4+0.2\left(\frac{T_w}{T_{pc}}-1\right),T_b\le T_{pc}\le T_w\\ 0.4+0.2\left(\frac{T_w}{T_{pc}}-1\right),1\le \frac{T_b}{T_{pc}}\le 1.2\end{array}\right)\end{array}$(4)Nub=0.0069Re0.9bPrb¯¯¯¯¯¯0.66(ρwρb)0.43(1+2.4Dx)$N{u}_b=0.0069R{e}_b^{0.9}{\overline{P{r}_b}}^{0.66}{\left(\frac{{\rho }_w}{{\rho }_b}\right)}^{0.43}\left(1+2.4\frac{D}{x}\right)$(5)Nub=2.0514Re0.928bPr0.742b(ρwρb)1.305(μwμb)−0.669(c¯pcp,b)0.888(q+)0.792$N{u}_b=2.0514R{e}_b^{0.928}P{r}_b^{0.742}{\left(\frac{{\rho }_w}{{\rho }_b}\right)}^{1.305}{\left(\frac{{\mu }_w}{{\mu }_b}\right)}^{-0.669}{\left(\frac{{\overline{c}}_p}{c_{p,b}}\right)}^{0.888}{\left(q^{+}\right)}^{0.792}$(6)Nub=0.0061Re0.904bPrb¯¯¯¯¯¯0.684(ρwρb)0.564(μwμb)−0.184$N{u}_b=0.0061R{e}_b^{0.904}{\overline{P{r}_b}}^{0.684}{\left(\frac{{\rho }_w}{{\rho }_b}\right)}^{0.564}{\left(\frac{{\mu }_w}{{\mu }_b}\right)}^{-0.184}$(7)

图7展示了各个关联式的验证结果，综合来看，Gupta、Kim和Jackson关联式与计算数据误差最大，大部分结果的验证误差均大于30%。Bishop关联式有部分点的验证误差大于30%，Mokry关联式在小热通量和大质量流速时与计算数据吻合较好，验证误差在30%以内，但其他工况与计算数据的符合较差，验证误差均大于30%。对于Fan关联式，除高质量流速（case 9）和大管径（case 10）的工况的部分数据外，绝大多数数据验证结果误差均在30%以内，预测精度较高。根据2.1节的分析，传热恶化的发生主要由近壁区黏度和密度的变化导致的，而Fan关联式引入了壁面密度和黏度的修正，因此该关联式精度较高，推荐使用Fan关联式来预测单侧加热方形通道内超临界水的换热性能。

图7

图7   关联式计算与数值模拟结果的比较

Fig.7   Comparison between the correlations and present numerical results

2.3 不同结构强化换热效果比较

为防止竖直向上流动工况出现传热恶化，本节对比了2.1节所述多种强化换热结构的强化换热效果。图8为G=500 kg·m-2·s-1时不同强化换热结构的加热底面平均温度随加热长度的变化。可以看出，对传热恶化工况，5种强化换热结构均可降低壁面温度，其中，四种凹槽结构和并排的双通道均可以避免传热恶化，显著降低壁温。并排双通道、非对称倒角下游凹槽、倒角凹槽、非对称倒角上游凹槽和对称凹槽，在光滑通道传热恶化点处的壁温下降幅度分别为18.55%、16.12%、14.76%、14.65%和13.29%，强化换热效果为：并排双通道>非对称倒角下游凹槽＞倒角凹槽>非对称倒角上游凹槽>对称凹槽，上下双通道强化换热效果最差，不能避免传热恶化的产生。

图8

图8   加热底面平均温度随加热长度的变化

Fig.8   Average temperature on the heated surface against heating length

为综合评价强化换热结构的换热能力和阻力损失，引入等泵功条件下综合换热评价因子PEC。

PEC=Nu/Nu0(f/f0)1/3$PEC=\frac{Nu/N{u}_0}{(f/f_0{)}^{1/3}}$(8)

图9、图10分别表示了竖直向上流动不同强化换热结构Nu/Nu0与PEC随质量流速的变化。可以看出非对称倒角下游凹槽综合换热性能最优，非对称倒角下游凹槽在质量流速G=500~1000 kg·m-2·s-1时，PEC有微小的下降，但维持在一个较高的水平，且Nu/Nu0比值维持较高水平，说明在该流量区间，随着流量增加非对称倒角下游凹槽对流体的扰动增加，相较于光滑管，强化传热效果能维持较高水平。对称凹槽与倒角凹槽结构也具有一定的综合强化换热性能，并排双通道仅仅在质量流速G=500 kg·m-2·s-1时有一定的强化效果，其他情况下PEC<1，综合换热性能较差。除了上下通道，其他通道在G=500 kg·m-2·s-1时PEC能达到最大数值，是因为在该流量下，强化结构对流体的扰动影响最大，Nu/Nu0达到最大值，此时强化传热效果最佳。上下双通道的Nu/Nu0较小，强化换热效果有限，其PEC在所有质量流速下均小于1，综合换热性能不如光管。

图9

图9   不同强化换热结构的Nu/Nu0

Fig.9   Nu/Nu0 of different enhanced structures

图10

图10   不同强化换热结构的PEC

Fig.10   PEC of different enhanced structures

综上所述，在q=1000 kW·m-2、G=200~1500 kg·m-2·s-1的工况下，双通道结构换热效果有限，不适用于超临界水的强化换热；凹槽结构均具有一定的强化换热效果，其中非对称倒角下游凹槽结构强化换热效果最佳。

2.4 强化换热机理分析

本节对凹槽结构强化换热的机理进行详细分析。图11为G=500 kg·m-2·s-1时四种凹槽结构附近切面的流线图，可以看到凹槽区域产生明显的旋涡，导致该区域的流动分离，从而破坏边界层的发展，增加流体的扰动。流体在凹槽下游位置与加热壁面再接触，使得该处边界层变薄，从而强化换热。由于旋涡区域仅出现在凹槽结构内部，对主流流动影响较小，因此不会带来较大的压降。值得注意的是，凹槽内部的旋涡虽然增强了流体的扰动，但由于其速度较小，因此换热能力较弱，导致凹槽内部区域换热效果较差，强化的部位只发生在凹槽下游的再附着点附近。因此，在保证流动分离的前提下，应尽可能减小旋涡强度以降低阻力。对称凹槽结构由于壁面结构的突变，产生了较强的流动分离，导致在凹槽区域出现了一大一小两个旋涡，如图11（a）所示。较大的旋涡可以使主流产生流动分离从而强化换热，但较小的旋涡对换热并没有积极作用，反而由于增强了流体的扰动而增大该处的流动阻力。与对称凹槽相比，倒角凹槽壁面附近流线更为顺滑，消除了凹槽区域较小的旋涡，使得流阻降低。同时，在凹槽下游部位，再附着区域由尖角变为平滑过渡的圆角，相当于增加了强化区域，因此倒角凹槽比对称凹槽具有更高的传热性能。非对称倒角上游凹槽结构上游陡峭而下游平缓，产生了一大一小两个旋涡，由于有利于强化传热的旋涡产生于下游，而下游的区域较为狭窄，所以非对称倒角上游强化传热效果一般。非对称倒角下游凹槽结构壁面下游陡峭而上游平缓，增强了流体的扰动，产生了一大三小的四个旋涡，大的旋涡对倒角凹槽结构影响范围更大，强化换热效果更为明显，小的旋涡则增加了流动阻力，综合来看大的旋涡影响范围更广，因此非对称倒角凹槽结构强化换热效果更好，推荐使用非对称倒角凹槽结构以改善传热恶化。

图11

图11   凹槽附近剖面流线图

Fig.11   Streamlines around various grooves

为进一步揭示传热强化机理，引入Guo等[32]提出的场协同原理。一方面是温度梯度场与流场的协同关系，即温度梯度与速度的协同角β，协同角β越小，对流传热能力越好；另一方面是速度矢量与速度梯度矢量的协同关系，即速度与速度梯度的协同角α，协同角α越大，流动阻力越小。

图12、图13均为沿着流动的中心面剖面图，分别表示G=500 kg·m-2·s-1时对称凹槽、倒角凹槽和非对称倒角凹槽在Z=0.25 m附近的协同角β与α的分布。可以看出，协同角β值较小的范围：非对称倒角下游凹槽>倒角凹槽>对称凹槽>非对称倒角上游凹槽，说明非对称倒角下游凹槽结构的传热能力最强。此外，对称凹槽结构中，β值较小的区域主要集中在凹槽内部，但由于该区域流速较小，换热能力较差，因此对主流换热能力的强化效果有限。而倒角凹槽与非对称凹槽结构中β值较小的区域主要集中在凹槽与主流区域交界处，该处的换热强化有利于凹槽区域和主流之间的热量交换，且非对称凹槽β值较小的区域向主流区沿伸幅度大于倒角凹槽，因此，非对称倒角凹槽结构更有利于传热的强化。

图12

图12   凹槽附近的协同角β云图

Fig.12   Distribution of synergy angle β around various grooves

图13

图13   凹槽附近的协同角α云图

Fig.13   Distribution of synergy angle α around various grooves

由图13（a）可以看出，对称凹槽上游和下游部位存在较大的低α区域，这是由于对称凹槽结构尖角的存在，流体在流动分离点和再附着点附近产生比较强烈的流动方向转变[图11（a）]，使得该处流动阻力严重增加。对凹槽结构采用倒角后，可以使流动更平滑，避免较强的流动分离及流动方向转变，从而显著改善这两个区域的协同性，使局部阻力显著降低。采用倒角的结构中低α区域影响区域：非对称倒角下游凹槽>倒角凹槽>非对称倒角上游凹槽，说明非对称倒角下游凹槽结构流动阻力最大，非对称倒角下游凹槽结构由于下游壁面变陡，增强了流体的扰动，强化换热的同时也增加了局部阻力，非对称倒角上游凹槽虽然降低了局部阻力但强化换热效果一般，综合PEC考虑，非对称倒角下游凹槽结构强化换热效果最佳。

3 结 论

建立了超临界水在单侧加热方通道内流动传热计算模型，基于数值模拟研究，分析了超临界流体常用换热关联式对超临界水在单侧加热方管中的适用性，计算结果可进一步为发动机再生冷却通道的设计提供参考。同时详细分析了传热恶化产生的机理，比较不同强化换热结构（双通道和凹槽）的优劣及机理分析。主要结论如下。

（1）引起单侧加热方管内超临界水流动传热恶化主要有以下两点原因：①黏性底层内黏度增大使得边界层局部增厚；②黏性底层与缓冲层内密度减小导致湍动能降低。上述耦合影响导致传热恶化的出现。

（2）综合比较了六种关联式，发现Fan关联式与计算结果符合程度最高。传热恶化的发生主要由近壁区黏度和密度的变化导致，而Fan关联式引入了壁面密度和黏度的修正，因此推荐采用Fan关联式预测单侧加热方形通道内超临界水的换热性能。

（3）在传热恶化产生的区域，综合对比了双通道和凹槽结构的强化传热效果，发现双通道结构强化传热效果有限，不适用于超临界水的强化传热；凹槽结构具有较好的强化传热效果，且引起的压力损失小，其中非对称倒角下游凹槽结构强化传热效果最佳。

（4）场协同原理分析发现，对于非对称倒角下游的凹槽结构，其β值较小的区域占比大，且低α区域增加不明显，说明其强化传热的同时，流动阻力增加不显著，因此非对称倒角下游凹槽结构综合传热性能最佳。

符 号 说 明