高温热辐射效应下相变渗透冷却的建模与数值研究

《International Journal of Heat and Mass Transfer》：Modeling and numerical study of phase change transpiration cooling under high temperature thermal radiation effects

　　

## 摘要

多孔结构中的相变渗流冷却是高热通量环境下一种有效的热防护方法，但在热辐射和非均匀加热条件下其瞬态传输行为仍缺乏充分理解。本文建立了物理模型，并基于OpenFOAM开发了数值求解器，研究多孔结构中液冷剂渗流冷却的瞬态传热传质特性。研究了热辐射、孔隙参数和非均匀热通量分布的影响。结果表明，热辐射是高温渗流冷却中一种重要的传热方式。吸收系数增大导致辐射穿透深度减小，辐射能集中在多孔结构的上部区域，从而导致冷却启动延迟、峰值温度升高以及内部温度梯度增大。多孔泡沫的每英寸孔数（PPI）对传热传质过程有显著影响。PPI增大促进两相区内的液体积聚，并加剧间歇性蒸汽排出行为。相比之下，低PPI泡沫结构有利于更平稳的蒸汽排出，并更有效地抑制瞬态温度峰值。在非均匀热通量作用下，高热通量侧的强烈蒸发引起横向压力梯度，驱动蒸汽向低热通量侧迁移，导致冷却能力与热载荷之间出现局部不匹配。这些结果加深了对多孔结构中瞬态相变渗流冷却的理解，为复杂热条件下渗流冷却系统的设计提供了有用的指导。

## 1. 引言

为解决高速飞行器在高热通量和长时间飞行条件下面临的热防护挑战，近年来主动热防护技术受到了越来越多的关注。在现有方案中，渗流冷却被广泛认为是一种有前景的解决方案，因为它可以通过冷却剂在多孔介质中的渗透、流体-固体传热以及相变吸热的共同作用来降低热载荷[1], [2]。然而，当前研究进展与实际工程应用之间仍存在较大差距。主要障碍之一在于多孔结构内部两相流、相变、共轭传热和热辐射传输之间的强耦合作用，底层的传热传质机理尚未被完全理解。因此，为推动渗流冷却技术的成熟，必须进一步阐明相关的传输过程，并开发可靠的物理模型和预测方法，用于结构优化和工程设计。

许多研究者使用气体冷却剂[3], [4]和液体冷却剂[5], [6], [7], [8], [9], [10], [11]进行了实验，验证了渗流冷却的可行性。Gülhan等人[12]在Ma 6风洞中比较了空气、氩气和氦气的冷却能力，发现在相同质量流量下氦气提供了最佳冷却效率。由于气液相变具有巨大的潜热，液体冷却剂相比气体冷却剂表现出更好的热防护性能。水和含水溶液是常用的液体冷却剂。在液体冷却剂的驱动方式方面，目前的主要方案为泵强制驱动和毛细力自驱动。

对于强制驱动的液体冷却剂供给方式，Wang和Zhao等人[13], [14]设计了一个半径为15°的楔形结构，并在风洞中进行了实验以验证其冷却能力。遗憾的是，高压区同时也是高热通量区，楔形表面的非均匀压力分布导致在强制驱动方式下大部分冷却剂流向低压区域[15]。这导致冷却剂流量与气动热通量之间出现不匹配。为克服这一问题，Wu和Wang等人[16]设计了一种分段孔隙率布局的渗流矩阵，在楔形前端使用渗透性更大的烧结316L合金多孔材料。这种非均匀渗透率布局缓解了冷却剂流量不匹配问题。Zhang和Huang等人[17]制备了微孔陶瓷材料作为渗流矩阵，并在风洞中进行了实验。他们发现实验过程中陶瓷材料的表面温度始终保持在360 K。Gulli等人[18]在弧风洞实验中采用碳-碳陶瓷材料作为渗流矩阵，验证了渗流冷却对受热表面的防护效果。Rocher等人[19]采用二硼化锆陶瓷作为渗流矩阵，在风洞中验证了渗流冷却防止表面氧化的效果。

清华大学蒋教授课题组提出了毛细力驱动的冷却剂供给方法，并对此进行了广泛研究。Huang和Jiang等人[20]采用烧结金属作为渗流矩阵，研究了烧结金属粒径对阻力和毛细力的影响，并提出了满足冷却剂自驱动流动的合适粒径范围。在此基础上，他们利用增材制造技术制备了三通道结构，用作冷却剂的自驱动流动结构。与烧结金属相比，该结构具有灵活可控的优势[21]。为减少气流对液滴的夹带并提高冷却剂利用率，他们在烧结金属外侧铺设了一层金属泡沫，冷却剂消耗量减少了约40%[22]。

在相变渗流冷却的数值研究中，多孔介质中的气液相变过程具有挑战性。目前可用的气液相变模型分为分离相模型（SPM）和两相混合模型（TPMM）。在分离相模型中，需要求解多孔介质内每种流体相的流动、传热和相变方程。虽然该方法能捕捉更多物理细节，但计算成本更高，且模型精度依赖于相间相互作用的建模。两相混合模型由Wang和Beckerman[23], [24]提出，因其计算成本低而被广泛应用。Su和He等人[15]采用TPMM结合商业软件Fluent，数值研究了高速流场与多孔介质中气液相变过程之间的耦合。由于原始TPMM中相界面附近对流系数的剧烈变化导致数值不稳定，Su和Wang等人[25]用混合温度格式能量方程替换了混合焓格式能量方程，从而改善了数值解的收敛性。Hu和Jiang等人[26]提出了一种新的动力学焓格式能量方程，有助于稳定求解过程并提高收敛性。

尽管TPMM及其改进版本已被广泛应用于解决多孔介质中的气液相变问题，但不可压缩流体假设使其不适用于气相密度变化显著的场景。Dong和Wang等人[27]提出半混合模型（SMM），用混合焓替换SPM中的能量方程，并研究了非均匀热通量下相变渗流冷却的特性。Wang等人[28]基于开源软件平台OpenFOAM中的SPM，开发了多孔介质中瞬态相变和流动的数值求解器，并研究了渗流矩阵中的气液相变过程。尽管SMM和SPM可以应用于更广泛的情况，但未考虑高温条件下的辐射传输现象。Zhang等人[29]利用Rosseland辐射传输模型报道了辐射传输对气体渗流冷却的影响。他们发现辐射传热使壁温升高，从而降低了整体冷却效率，并增强了从高温气体到壁面的总传热。因此，有必要考虑辐射传输对渗流冷却的影响。

本文基于OpenFOAM-v2406开发了一个考虑多孔介质中流动、辐射传输和相变的数值求解器框架。基于该求解器，对高热通量作用下的液相变渗流冷却进行了数值研究。本文的组织结构如下：第2节介绍简化的主动热防护结构及其涉及的物理过程。第3节建立物理模型，第4节给出数值实现和模型验证。第5节讨论液相变渗流冷却的特性，第6节得出结论。

## 2. 几何模型与问题描述

考虑一种双层多孔热防护结构，如图1所示。采用双层配置是为了代表高热通量环境下实际渗流冷却热防护系统中常用的功能分层设计。结构A是一种耐高温烧结多孔材料，孔隙率为?=0.3，作为直接暴露于恶劣热环境的热防护层。该层设计用于承受高热载荷并允许蒸汽在多孔基质内充分进行热交换，从而在排出前达到过热状态。结构B是一种高孔隙率多孔泡沫，?=0.9，作为冷却剂传输层以降低流动阻力、促进蒸汽排出并维持稳定的冷却剂供给。A和B的厚度分别为3 mm和7 mm，宽度为20 mm。结构A的上表面被加热，冷却剂入口位于结构B的底部。结构B底部延伸0.5 mm以定义冷却剂入口边界。


在上述问题中，涉及多孔介质中的传热、流动和相变，由于固体结构温度高，还考虑了辐射传热。因此，数值求解该问题具有挑战性，我们将在下一节讨论上述物理过程的建模。

## 3. 物理模型与方法

多孔介质中流动和传热的物理模型可分为连续尺度模型和孔隙尺度模型。其中，连续尺度模型因计算成本低而被广泛采用。本节在连续尺度上对多孔介质内的流动、传热和相变进行建模。建模过程中采用以下假设：

(1) 两相区内气相和液相的速度相同；
(2) 气液界面处于局部热力学平衡状态；
(3) 流体与固相之间存在有限温差，在两相区内蒸汽和液体处于热平衡；
(4) 冷却剂为单一组分流体。

关于第一个假设，气液相速度相同是采用均匀混合流近似的本连续尺度模型。求解单一混合动量方程，不显式考虑气相和液相之间的相对滑移速度。该假设仅适用于液体和蒸汽共存的两相区域。在高温干区，液相体积分数变为零，动量方程描述的是蒸汽传输而非两相同时运动。类似的处理方法已被采用于多孔wick中液-汽相变、多孔通道干燥过程以及金属泡沫沸腾传热的先前研究中[30], [31], [32]。然而，局部液-汽滑移的可能影响将在未来工作中使用包含相间滑移关系的混合动量模型或两流体模型进一步研究。此外，辐射传输模型的假设在第3.3节单独讨论。

### 3.1. 相分数方程

相分数方程描述了多孔介质内各相的传输过程。液相和汽相的相分数方程为[33]：

(1) ? ?(ρ?α?)/?t + ?·(?ρ?α?u) = ??

其中下标i为相指数，表示液相或汽相，?为孔隙率，α?和ρ?分别为相分数和密度，u为固有平均速度，??为相变速率。考虑质量守恒，汽相的质量增加应等于液相的质量减少。

### 3.2.

# 动量方程与能量守恒方程

基于第一个假设，我们可以推导出多孔介质内流体相的动量方程，即[33]，

(2) $\frac{\partial(\varphi \rho_m \mathbf{u})}{\partial t} + \nabla \cdot (\varphi \rho_m \mathbf{u}\mathbf{u}) = -\varphi \nabla p + \alpha_l \varphi \nabla p_c + \nabla \cdot (\varphi \boldsymbol{\tau}) + \varphi \rho_m \mathbf{g} - \varphi^2 \sum_i \left(\frac{\alpha_i \mu_i}{K K_{r,i}}\right) \mathbf{u}$

其中 $\rho_m$ 为相分数加权平均密度，$p$ 和 $p_c$ 分别为流体相的热力学压力和毛细管压力，下标 $l$ 代表液相，$\boldsymbol{\tau}$ 为流体相的粘性应力，$\mu_i$ 为相 $i$ 的动力粘度，$\mathbf{g}$ 为重力加速度，$K$ 为多孔介质的固有渗透率，$K_{r,i}$ 为相 $i$ 的相对渗透率。毛细管压力和相对渗透率的建模将在第3.4节中介绍。

在连续尺度上，多孔介质中流体与固体之间的传热模型可以分为局部热非平衡（LTNE）模型和局部热平衡（LTE）模型。在多孔介质中，由于传热发生在液相、气相和多孔固体相之间，这一问题变得更加复杂。根据第三个假设，流体相的能量方程为[33]：

(3) $\left(\varphi \sum_i \alpha_i \rho_i C_{p,i}\right) \left(\frac{\partial T_f}{\partial t} + \mathbf{u} \cdot \nabla T_f\right) = \nabla \cdot \left(\sum_i \alpha_i \varphi k_i \nabla T_f\right) + \dot{m}_l h_{lg} + S_{r,l} + S_{sf}$

其中 $T_f$ 为流体温度（液相和气相），$C_{p,i}$ 为相 $i$ 的定压比热容，$k_i$ 为相 $i$ 的热导率，$h_{lg}$ 为潜热，$S_{r,l}$ 为作用于液相的辐射源项，$S_{sf}$ 为流体相与多孔固体相之间对流换热的源项，其计算公式为：

(4) $S_{sf} = \sum_i (\alpha_i h_{s,i}) A_{sf} (T_s - T_f)$

其中 $h_{s,i}$ 为多孔固体相与相 $i$ 之间的传热系数，$A_{sf}$ 为多孔介质的比表面积，$T_s$ 为多孔固体相的温度。在连续尺度上，必须补充传热系数以保证传热方程的完整性，这将在第3.4节中讨论。

对于多孔介质中的固体相，能量方程为：

(5) $(1-\varphi) \rho_s C_{p,s} \frac{\partial T_s}{\partial t} = \nabla \cdot \left((1-\varphi) k_s \nabla T_s\right) + S_{r,s} - S_{sf}$

其中 $\rho_s$、$C_{p,s}$、$k_s$ 和 $T_s$ 分别为多孔固体相的密度、比热容、热导率和温度，$S_{r,s}$ 为作用于多孔基质的辐射源项。

## 3.3 辐射传输方程

辐射传热作为传热方式之一，在传热过程中起着重要作用，尤其是在高温情况下。在单相流中，气体被视为不参与辐射传输过程的透明介质。在连续尺度上，多孔结构内的热辐射传输可以视为等效半透明介质内的辐射传输[34]、[35]。因此，建立同时包含液相和多孔基质的连续尺度辐射传输模型需要以下假设：

(1) 气相不参与辐射传热。

(2) 不考虑散射效应。

(3) 液相被视为半透明介质。

(4) 液相和固体相的辐射吸收系数彼此独立。

(5) 固体相为灰体。

这些假设用于构建当前多相多孔介质的等效连续尺度辐射模型。气相被视为不参与，因为多孔结构的特征尺寸仅为厘米量级，在当前条件下孔内空间中的气体为光学薄层。因此，其辐射吸收预计远弱于多孔固体基质和液态水的辐射吸收。忽略散射效应是因为散射主要重新分配辐射方向和空间分布，而吸收决定辐射能量向热能的转换。此外，在包含固体骨架、液态水和水蒸气的非均质多孔结构中，散射系数和相函数难以准确确定。因此，对于当前连续尺度模型，忽略气相辐射和散射是合理的。固体相的灰体假设也作为多孔固体骨架辐射特性的等效处理。已有研究表明，对于开孔泡沫状多孔材料，其辐射特性与孔隙率、孔径和骨架形态密切相关，而波长依赖性相对较弱[36]、[37]、[38]。对于烧结多孔材料，不透明多孔陶瓷层的辐射特性也被报道为随波长变化较小，支持使用灰体近似[39]。因此，在当前连续尺度辐射模型中，将固体相视为灰体被认为是合理的。

基于这些假设，一阶球谐近似（P1）下的辐射传输方程为[40]：

(6a) $\nabla \cdot \left(\Gamma_\lambda \nabla G_\lambda\right) + 4\pi \left[\alpha_l \kappa_{l,\lambda} n_{l,\lambda}^2 I_{b,\lambda}(T_f) + \kappa_s n_s^2 I_{b,\lambda}(T_s)\right] = (\alpha_l \kappa_{l,\lambda} + \kappa_s) G_\lambda$

(6b) $\Gamma_\lambda = \frac{1}{3}(\alpha_l \kappa_{l,\lambda} + \kappa_s)$

(7a) $S_{r,l} = \alpha_l \int_0^\infty \kappa_{l,\lambda} \left[G_\lambda - 4\pi n_{l,\lambda}^2 I_{b,\lambda}(T_f)\right] d\lambda$

(7b) $S_{r,s} = \kappa_s \int_0^\infty \left[G_\lambda - 4\pi n_s^2 I_{b,\lambda}(T_s)\right] d\lambda$

其中 $G_\lambda$ 为光谱入射辐射，$I_{b,\lambda}$ 为光谱黑体辐射强度，$\kappa_{l,\lambda}$ 和 $n_{l,\lambda}$ 分别为液相的等效光谱辐射吸收系数和光谱折射率，$\kappa_s$ 和 $n_s$ 分别为多孔材料固体相的等效辐射吸收系数和折射率。

## 3.4 模型的封闭

在连续尺度上，需要额外的模型来封闭流体的控制方程。在式(2)中，各相的相对渗透率 $K_{r,i}$ 采用Brooks-Corey模型[41]计算：

(8) $K_{r,i} = \left(\frac{\max(\alpha_i, \alpha_{res})}{\alpha_{res}}\right)^{3/2}$

其中 $\alpha_{res}$ 为残余相分数，取值为0.001。毛细管压力 $p_c$ 由以下关联式[42]计算：

(9a) $p_c = \frac{\sigma_l \cos(\theta_c)}{(\varphi/K)^{0.5}} J(\alpha_v)$

(9b) $J(\alpha_v) = 1.417 \alpha_v - 2.12 \alpha_v^2 + 1.263 \alpha_v^3$

其中 $\sigma_l$ 为液相的表面张力系数，$\theta_c$ 为液相与固体相之间的接触角，取值为60°，$J$ 为关联式[24]。

式(4)中的传热系数根据多孔结构的特性进行评估。对于烧结多孔结构A，采用Wakao等人提出的关联式[43]：

(10a) $h_{s,i} = \left(2 + 1.1 Re_i^{0.6} Pr_i^{0.33}\right) \frac{k_i}{d_c}$

其中 $d_c$ 为特征长度，等于烧结颗粒的直径，$Re_i$ 和 $Pr_i$ 分别为相 $i$ 的雷诺数和普朗特数。对于多孔泡沫，广泛使用以下关联式[44]：

(10b) $h_{s,i} = \left(32.504\varphi^{0.38} - 109.94\varphi^{1.38} + 166.65\varphi^{2.38} + 86.98\varphi^{3.38}\right) \frac{k_i}{A_{sf} d_c^2 Re_i^{0.438}}$

其中 $d_c$ 等于泡沫的孔径。在双相区，多孔固体相与气相之间的传热系数由式(10b)评估，而多孔固体相与液相之间的传热系数由Rohsenow提出的沸腾换热关联式[45]确定：

(11) $h_{s,i} = \mu_l \left(\frac{C_{p,l}}{C_{wl} h_{lg} Pr_l^{1.7}}\right)^{1/3} \frac{(T_s - T_{sat})^2}{h_{lg}^2} \sqrt{\frac{g(\rho_l - \rho_v)}{\sigma_l}}$

其中 $T_{sat}$ 为饱和温度，$C_{wl}$ 为经验常数，取值为0.006[28]。

此外，需要相变公式来计算相变速率 $\dot{m}_l$。基于动力学理论的Schrage模型提供了液-汽界面蒸发和凝结的理论描述[46]。然而，其应用需要指定适应系数，目前仍缺乏可靠且通用的确定准则，尤其是对于复杂多孔介质中的相变。因此，在当前连续尺度多孔介质模型中，采用基于能量平衡的质量传输速率公式。根据第二个假设，气-液界面处局部热力学平衡成立，双相区内的流体温度等于局部饱和温度。在此条件下，扣除显热变化、对流传输、热传导、辐射吸收和固-流体换热后，供给双相区的能量转化为液-汽相变的潜热。因此，可以在双相区内复用式(3)来计算相变速率，如下所示：

(12) $\dot{m}_l = \frac{1}{h_{lg}} \left[\left(\varphi \sum_i \alpha_i \rho_i C_{p,i}\right) \left(\frac{\partial T_f}{\partial t} + \mathbf{u} \cdot \nabla T_f\right) - S_{r,l} - S_{sf} - \nabla \cdot \left(\varphi \sum_i \alpha_i k_i \nabla T_f\right)\right]_{T_f = T_{sat}}$

与焓方法相比，当前公式通过局部能量平衡显式确定液-汽质量传输速率，而不是通过等效焓或热容隐式处理相变。与Lee型模型不同，它不需要经验相变松弛系数。在数值实现中，双相区内的流体温度被约束在局部饱和温度，质量传输速率由离散化的流体能量方程获得，通过将剩余能量残差转化为潜热来实现。

## 3.5 边界条件

如图1所示，结构A的上边界为开放边界，蒸汽由此流出，固体相上施加热通量，并考虑向周围环境的辐射传热。结构B的下边界为冷却剂入口。结构A和B的侧边界视为对称。对于辐射传热，所有壁面假设为漫发射灰体表面。上述边界条件的数学表达式如下：

对于结构A的上边界，边界条件为：

(13) $\frac{\partial u}{\partial n} = 0,\ \frac{\partial \alpha_l}{\partial n} = 0,\ p = p_0,\ \frac{\partial T_f}{\partial n} = 0,\ (1-\varphi)k_s \frac{\partial T_s}{\partial n} + q_{r,w} = q_0 - \varepsilon \sigma (T_s^4 - T_0^4)$

其中 $p_0$ 和 $T_0$ 分别为周围静压力和温度，$q_{r,w}$ 为壁面处的辐射热通量，$q_0$ 为作用于固体相的热通量，$\varepsilon$ 为壁面的发射率，取值为0.8。值得注意的是，对于多孔固体温度边界条件，左侧的导热热通量和内部P1辐射通量代表传输到多孔结构中的能量，而右侧的项代表外部施加的热通量和壁面向环境的辐射热损失。冷却剂入口处的边界条件为：

(14a) $\rho_l |\mathbf{u} \cdot \mathbf{n}| = q_{m,in},\ \alpha_l = 1$

(14b) $\frac{\partial p}{\partial n} = -\frac{1}{\varphi} \mathbf{n} \cdot \left[\frac{\partial(\varphi \rho_m \mathbf{u})}{\partial t} + \nabla \cdot (\varphi \rho_m \mathbf{u}\mathbf{u}) - \nabla \cdot (\varphi \boldsymbol{\tau}) - \varphi \rho_m \mathbf{g} + \varphi^2 \sum_i \left(\frac{\alpha_i \mu_i}{K K_{r,i}}\right) \mathbf{u}\right] - \alpha_l \frac{dp_c}{d\alpha_v} \frac{\partial \alpha_l}{\partial n}$

(14c) $T_f = T_{in},\ \frac{\partial T_s}{\partial n} = 0$

其中 $\mathbf{n}$ 为边界的外法向单位矢量（指向域外），$q_{m,in}$ 为冷却剂的质量通量，$T_{in}$ 为冷却剂的入口温度。结构A和B的侧边界设为对称。对于辐射，结构A的上边界和冷却剂入口的条件为：

(15) $-\Gamma_\lambda \frac{\partial G_\lambda}{\partial n} = \varepsilon(4-2\varepsilon)\left[G_\lambda - 4\pi n_i^2 I_{b,\lambda}(T_s)\right],\ q_{r,w} = -\int_0^\infty \Gamma_\lambda \frac{\partial G_\lambda}{\partial n} d\lambda$

其中 $T_w$ 为壁面温度，$n_i$ 为相邻相的折射率，其值取决于边界相邻的相。

## 4. 数值策略与验证

### 4.1 数值实现

气相、液相和多孔固体相的能量方程采用有限体积法离散和求解。光谱P1方程采用光谱带方法数值求解，即将连续光谱域离散为有限数量的波段。在本节中，我们讨论相分数方程和压力方程数值求解的细节。

对于相分数方程(1)，将其改写为以下形式：

(16) $\varphi \frac{\partial \alpha_i}{\partial t} + \nabla \cdot (\alpha_i \varphi \mathbf{u}) = -\frac{\alpha_i}{\rho_i} \varphi \frac{D\rho_i}{Dt} + \frac{\dot{m}_i}{\rho_i}$

将式(16)求和，并注意到 $\sum_i \alpha_i \equiv 1$，得到：

(17) $\nabla \cdot (\varphi \mathbf{u}) = -\sum_i \frac{\alpha_i}{\rho_i} \varphi \frac{D\rho_i}{Dt} + \left(\frac{1}{\rho_l} - \frac{1}{\rho_v}\right) \dot{m}_l$

在式(16)右侧添加和减去 $\alpha_l \nabla \cdot (\varphi \mathbf{u})$，利用式(17)，经过代数运算后，得到液相的守恒方程：

(18) $\varphi \frac{\partial \alpha_l}{\partial t} + \nabla \cdot (\alpha_l \varphi \mathbf{u}) = \alpha_l \alpha_v \left(\frac{1}{\rho_v} \frac{D\rho_v}{Dt} - \frac{1}{\rho_l} \frac{D\rho_l}{Dt}\right) + \left[\frac{1}{\rho_l} + \alpha_l \left(\frac{1}{\rho_v} - \frac{1}{\rho_l}\right)\right] \dot{m}_l + \alpha_l \nabla \cdot (\varphi \mathbf{u})$

在推导压力守恒方程之前，引入参考压力 $p_{rgh} = p - \rho_m \mathbf{g} \cdot \mathbf{r}$，我们得到以下关系：

(19) $\nabla p = \nabla p_{rgh} + \rho_m \mathbf{g} + \mathbf{g} \cdot \mathbf{r} \nabla \rho_m$

其中 $\mathbf{r}$ 为空间位置矢量。将式(19)代入式(2)并用有限体积法离散，得到：

(20a) $\mathbf{u}_P = \mathbf{u}_P^* - \frac{\varphi}{a_P} \nabla p_{rgh}$

(20b) $\mathbf{u}_P^* = -\sum_N \frac{a_N}{a_P} \mathbf{u}_N + \frac{b}{a_P} - \frac{\varphi}{a_P} \left(\mathbf{g} \cdot \mathbf{r} \sum_{f^n} \rho_{m,f^n} \mathbf{S}_{f^n} + \alpha_l \frac{dp_c}{d\alpha_v} \sum_{f^n} \alpha_{l,f^n} \mathbf{S}_{f^n}\right)$

其中 $a$ 为矩阵对角线上的系数，$b$ 为线性方程组的右端项，下标 $P$ 和 $N$ 分别表示当前单元和其相邻单元，$\mathbf{u}_P^*$ 为预测速度，$\mathbf{S}_{f^n}$ 为面积法向矢量，下标 $f^n$ 表示面心处的值。将式(20a)代入式(17)，最终得到以下压力方程：

(21) $\nabla \cdot (\varphi \mathbf{u}_P^*) - \nabla \cdot \left(\frac{\varphi^2}{a_P} \nabla p_{rgh}\right) = \dot{m}_l \left(\frac{1}{\rho_l} - \frac{1}{\rho_v}\right) - \sum_i \frac{\alpha_i}{\rho_i} \varphi \frac{D\rho_i}{Dt}$

考虑到传热和传质过程的耦合，以及辐射传输在多孔介质中带来的数值求解额外困难，我们引入内部迭代来解决这一问题。数值求解过程如图2所示。

*图2. 数值算法流程图。*

### 4.2 模型验证

在本节中，采用一维移动边界问题来验证模型的准确性。初始状态下，静止的蒸汽和水被相界面分隔。气相由过热壁面加热，液相侧的壁面温度固定在饱和温度，如图3所示。

*图3. Stefan问题的示意图。*

由于过热壁面的作用，液体持续蒸发，界面位置 $\delta_0(t)$ 和气相的无量纲温度场 $\Theta(y, t)$ 具有以下精确解：

(20) $\delta_0(t) = 2\eta \sqrt{\frac{k_v}{\rho_v C_{p,v}} t}$

(21) $\Theta(y, t) = \frac{T(y,t) - T_w}{T_{sat} - T_w} = \frac{1}{\text{erf}(\eta)} \text{erf}\left(\frac{y}{\delta_0 \eta}\right)$

其中 $\eta$ 为无量纲因子，由式(22)确定：

(22) $\eta \exp(\eta^2) \text{erf}(\eta) = \frac{C_{p,v}(T_w - T_{sat})}{\sqrt{\pi} h_{lg}}$

在此情况下，水和蒸汽分别选为液相和气相。

数值计算在0.1 MPa的饱和压力下进行，热壁面的过热度设定为10 K。计算网格的空间分辨率为5 μm。将水和蒸汽在饱和压力下的热物性参数代入式(22)，得到η = 0.0677745。图4展示了界面位置δ0以及y = 1 mm处的瞬态无量纲温度与数值结果的对比，表明当前的传质模型能够捕捉界面位置和温度的变化规律。

接下来，我们将考虑多孔介质内的沸腾传热问题。如图5所示，一个直径为40 mm、长度为100 mm的烧结青铜多孔介质，其底部由0.21 MW/m2的固定热通量加热。多孔介质的孔隙率和颗粒直径分别为0.318和0.2 mm。在多孔介质底部，水通过长度为10 mm的源项注入，源项的孔隙率设定为0.999。在初始时刻，水填充多孔介质的一半，温度设定为300 K。稳态下的固体温度数值结果与文献[26]中的实验数据如图6所示。可以看出，数值结果的温度与实验数据吻合，最大相对误差为4.87%。

为了验证辐射模型及其数值实现，引入了一个二维纯辐射传热基准问题。该基准问题由一个带有内管的半圆形腔体组成，其中填充有灰色、吸收-发射且不散射的介质，如图7所示。介质温度固定在T = 1000 K，所有固体壁面均作为黑体处理，Tw = 0 K。外半圆形壁面和内管的半径分别为1 m和0.2 m。考虑三种吸收系数κ = 1、3和10 m?1，以检验辐射模型在不同光学厚度下的性能。

图8比较了P1模型和DOM模型预测的外半圆形底部壁面上的归一化壁面辐射热通量分布与精确解[47]的结果。如图所示，P1和DOM模型均再现了壁面辐射热通量的对称分布，并捕捉到了吸收系数引起的变化。随着吸收系数的增大，辐射热通量分布变得更加均匀，表明参与介质的体积吸收作用更强。P1结果在当前光学厚度范围内与精确解表现出合理的吻合度，尽管由于P1模型是一种扩散型近似，无法完全解析边界附近强方向性的辐射，因此可以观察到局部偏差。DOM结果也与基准解吻合良好，为数值辐射处理提供了有用的交叉验证。此验证确认了本代码中使用的辐射传热方程和辐射边界处理是可靠的。

5. 结果与讨论

本节基于建立的模型，对图1所示的双层主动热防护结构内的传热和相变过程进行数值分析。计算网格总数为9000，最大网格尺寸为0.2 mm。该网格数量满足网格独立性要求。结构A和B的热物性参数选取自文献[48]、[49]，导热系数和比热容如图9所示。由于水安全且易于获取，选择水作为冷却剂。蒸汽被视为理想气体，蒸汽压方程由Antoine关联式给出。

泡沫的多孔结构如图10所示。在结构参数中，孔径dc和孔隙直径dp在决定泡沫整体性能方面起主导作用。在实际应用中，具有不同孔隙结构的泡沫材料通常以孔隙率和每英寸孔数（PPI）来表征。因此，明确孔隙参数对液体冷却剂传热特性和相变行为的影响至关重要。对于结构B，选择了三种具有代表性的孔隙构型进行研究。结构A和B的孔隙参数和等效辐射特性列于表1。值得注意的是，在干限条件下αl = 0时，PPI值为20、40和60的结构B的光学厚度分别为2.47、3.49和5.48。此外，当考虑液相时，液态水的强吸收可以显著增加有效光学厚度。例如，在液相比例为0.01时，有效光学厚度增加到约10.87–18.08，而在完全饱和时可达到842.47–1265.48。因此，结构B可被视为中等到大光学厚度的介质，在本连续介质尺度辐射传热计算中使用P1模型是合适的。泡沫结构比表面积和渗透率采用Dietrich等人[50]提出的关联式计算，等效吸收系数根据Li和Xia等人[51]建立的关联式确定。水在各光谱带中的吸收系数和折射率列于表2。

为了验证网格独立性，在热通量为1 MW/m2、冷却剂质量通量为150 g/(s·m2)的条件下，使用三种网格进行了数值模拟。Grid 1至Grid 3的网格单元数分别为6400、9000和12,000。选取t = 100 s时沿中心线的温度分布作为评价标准。如图11所示，三种网格产生的温度分布几乎相同，加热面温度分别为1406.75、1408.46和1408.65 K。Grid 2与最细网格之间的温度差仅为0.19 K。因此，选择Grid 2用于后续模拟。

时间步长根据Courant数约束动态调整。为评估数值结果对时间步长大小的敏感性，在Grid 2上进行了附加模拟，Comax = 0.4、0.6和0.9，其余条件保持不变。如图12所示，减小时间步长对平均加热面温度影响较小，t = 100 s时的温度分别为1408.22、1408.64和1409.25 K。三种情况之间的最大差异仅为1.03 K，对应的相对偏差约为0.073%。这些结果表明，在检查的Courant数范围内，数值解对自适应时间步长的进一步细化不敏感。因此，后续模拟中采用Comax = 0.9，以平衡数值精度和计算效率。

5.1. 辐射传热的影响

在多孔结构中气液相变的数值研究中，很少有研究人员考虑多孔介质相变过程中的辐射传热。因此，有必要进一步研究辐射传热对多孔结构中液体冷却剂传热和相变的影响。在本节中，结构A的表面热通量固定为q0 = 1 MW/m2，冷却剂质量通量设定为qm = 150 g/(s·m2)。后续分析针对PPI = 40的孔隙参数的结构B进行。

100 s时的液相分布和冷却剂温度场如图13所示。从图13(a)和(c)可以观察到，结构B中存在一个充满液体的区域，这里称为液池区域。当忽略热辐射效应时，液池区域位于比考虑热辐射时更靠近加热面的位置。然而，这并不意味着对加热面的冷却效应更强。由于热辐射是一种重要的传热模式，特别是在高温下，忽略热辐射效应会减少从加热面传递到结构中的热量。因此，液体冷却剂必须到达更靠近加热面的位置才能发生相变。图13(b)和(d)所示的冷却剂温度分布表明，结构B中液池区域上方存在一个两相区域，此处传热强度最高。随后，相变产生的蒸汽冷却剂在结构A中进行对流换热。被加热到过热状态后，蒸汽通过结构B表面的微孔排出。值得注意的是，虽然在加热面施加了均匀热通量边界条件，但液相和温度分布并不严格对称。这主要是由于多孔介质中液-汽相变、两相流和传热之间的非线性耦合所致。在蒸发/沸腾过程中，局部蒸汽生成改变了相分数、混合物密度、相对渗透率和流动阻力，从而影响液体补充和蒸汽排出路径。特别是动量方程中的对流项是非线性的，因此微小的局部速度扰动可能被放大并发展为优选蒸汽排出通道和不均匀蒸发前沿。因此，图13中观察到的不规则分布被认为是与蒸发/沸腾引起的两相流不稳定性相关的物理特征，而非数值伪影。在均匀加热条件下的相变渗透冷却过程中，也实验观察到了类似的表面温度波动和周期性加热/冷却行为[53]，这被归因于沸腾不稳定性以及蒸汽生成和液体补充之间的间歇性相互作用[9]。为确保数值收敛，在每个时间步内执行了额外的内循环迭代，所有控制方程的残差降至10??以下。

加热面的瞬时平均温度如图14(a)所示。温度分布表明，初始阶段表面温度迅速升高。随后，结构B中冷却剂相变产生的蒸汽为结构A提供对流冷却，导致表面温度逐渐下降并最终趋于稳定状态。由于热辐射增强了从结构A到结构B的传热，当考虑辐射传热时，结构A的加热面温度在达到稳态之前保持较低水平，最大温度差为68 K，发生在17.64 s。此外，从图14(a)可以观察到，在10 s之前热辐射效应对加热面温度的影响有限。由于烧结多孔材料的高吸收系数，辐射传热的贡献相对较小。根据基于Rosseland扩散理论的近似估算，1400 K时的辐射传热等效导热系数约为0.23 W/(m·K)，比结构A的本征导热系数低一个数量级。

结构A加热面的瞬态参数

如图14(b)所示，在约10 s时，加热表面的蒸汽质量流率呈现先快速下降后上升的阶段特征。这种行为可归因于毛细作用向上驱动冷却剂先发生相变，而产生的蒸汽引起的对流传热与热传导产生竞争。结果导致相变速率降低，随着传热过程的推进，相变速率又逐渐增大。此外，当忽略辐射传热时，蒸汽质量流率的快速下降和随后上升的起始时刻延迟了1.25 s，且蒸汽质量流率的下降更为明显。这是因为热辐射促进了热量向结构B的更快传递。同时，结构B较低的吸收系数使部分辐射能量能够穿透多孔泡沫层，直接作用于两相区域内的液相，从而提高液相温度。总体而言，辐射效应加速了结构B中冷却剂的相变过程，并使结构A的峰值温度降低。

图15展示了10 s和100 s时刻沿y方向中心线上的结构温度分布。考虑热辐射时，结构A的温度更低。一个显著特征是温度分布曲线在某一特定位置出现交叉。在该位置以上，考虑热辐射时温度更低；而在该位置以下，忽略辐射时温度反而更低。这一现象可由热辐射增强的传热效应来解释。相反，在交叉点以下区域，来自上部区域的辐射对结构产生净加热效应。

图15. 不同时刻y方向中心线上的结构温度

**5.2 孔结构参数的影响**

结构B的孔结构参数影响等效热物性、辐射特性及流动阻力，这些参数共同影响结构B内液态冷却剂的相变过程。图16展示了100 s时结构B内入射辐射能量的分布情况。结构A的上表面由恒定热流加热，辐射能量沿路径传播并因沿途吸收而衰减。根据表1，当结构B的PPI从20增加到60时，吸收系数从352.2 m?1上升至782.6 m?1。吸收系数的增加表明结构B内的辐射吸收能力增强，导致辐射能量衰减更为显著，其穿透效应随之减弱。可以预期，入射辐射能量的分布变化将显著影响结构内液态冷却剂的相变过程。

图16. 不同孔结构参数下结构B内入射辐射的分布

图17展示了不同孔结构参数下100 s时刻液相份额分布和冷却剂温度分布。如图17(a)、(c)和(e)所示，液池区域的厚度随PPI的增大而增加，且液池厚度与PPI之间并非线性关系。在PPI=20和40时，液池厚度保持在约2 mm左右；而在PPI=60时，液池厚度达到3 mm。液池厚度随PPI的变化与入射辐射能量的变化密切相关。如图16所示，PPI的增大缩短了入射辐射能量在结构B中的穿透距离，这种分布变化降低了同一位置的入射辐射能量，从而削弱了传热能力。液态冷却剂需要在更接近加热面的位置发生相变，导致液池厚度增加。此外，随着PPI增大，多孔泡沫的孔径减小，两相区域内的液相份额和冷却剂温度在不同PPI值下表现出不同的特征。如图17所示，在PPI=20和40时，两相区域呈现"羽毛状/指状"向上渗透的特征，液相分布和温度界面极不规则。当PPI增至60时，两相区域内的液相份额和冷却剂温度变得均匀。液相份额和冷却剂温度的变化与泡沫的流动阻力和传热能力的变化相关。

图17. 100 s时刻液相份额和冷却剂温度的分布

图18(a)给出了加热面蒸汽质量流率的瞬态变化。结构B加热面的蒸汽排出呈现明显的波动特征。PPI较高时，蒸汽排出表现出更强的间歇性，波动更为显著，且间隔相对较长。相比之下，较低的PPI则表现为波动周期较短、幅度较小，对应更频繁且相对均匀的蒸汽释放。

图18. 加热面的瞬态参数

图18(b)绘制了不同PPI值下结构A加热面的温度响应。表面温度响应可分为三个明显阶段。第一阶段为0至5 s，对应快速加热过程，此时加热面温度主要受结构A的热物性控制，而结构B的孔结构参数影响甚微。第二阶段为5至20 s，对应缓慢加热过程，表面温度逐渐达到峰值。随着传热推进，结构B温度升高，冷却剂在热传导和热辐射的共同作用下发生相变。在此阶段，结构B的孔结构参数对加热面温度有显著影响。峰值温度随PPI增大而升高，在PPI分别为20、40和60时，峰值温度分别为1413.02 K、1432.32 K和1499.75 K。第三阶段为20至100 s，对应准稳态阶段。由于蒸汽持续冷却作用，加热面温度整体呈下降趋势，并伴随明显波动。这种波动与液态补充、蒸汽生成和蒸汽排出之间的动态耦合有关。局部蒸发的增强促进了蒸汽的生成和排出，同时暂时降低了局部液体含量和补充量。当液态冷却剂重新进入相变区域时，蒸发再次加剧，导致反复的热振荡。此前关于相变透平冷却的实验研究中也报道了类似的非稳态温度波动和冷却不稳定性现象。Luan等人[9]观察到相变传热过程中明显的非稳态响应，He等人[10]将波动性热行为归因于蒸发/沸腾、液体补充和蒸汽排出之间的耦合效应。

图18(c)展示了固体与蒸汽之间的温差。随着蒸汽质量生成量的增加，局部热非平衡效应更加明显。然而，出口处固液温差整体上仍较小，在100 s的整个过程中最大值始终低于4.3 K。

图19展示了10 s和100 s时刻沿y方向中心线上的固体温度和冷却剂温度分布。不同时刻的温度分布进一步揭示了多孔结构参数对流体与固体基体之间共轭传热机制的影响。如图19(a)所示，在早期瞬态阶段（10 s），固体温度分布仅受PPI的微弱影响。同时，图18(a)表明此时相变产生的蒸汽质量流率仅为稳态值的约一半，因此结构内部的热响应主要受热物性控制。相比之下，图19(b)中冷却剂温度分布已开始表现出对PPI的敏感性，表明不同的多孔结构在冷却启动阶段对传热行为有显著影响。当系统在100 s逐渐接近准稳态时，多孔结构参数对传热的影响更加显著。对于PPI较高的泡沫结构，吸收系数的增大和辐射穿透深度的减小导致更多辐射能量沉积在结构B的上部区域。如图19(c)所示，结构B的温度随PPI增大而整体降低。同时，如图19(d)所示，较高PPI下结构B内两相区域的厚度明显减小。

图19. y方向中心线上的固体和冷却剂温度分布

为进一步揭示孔结构参数影响透平冷却性能的内在机理，本节对孔结构相关参数进行了解耦分析。分别独立改变渗透率、比表面积和有效吸收系数，同时保持其余两个参数不变，分析这些参数对峰值温度、冷却启动时间、液池厚度和蒸汽质量流率波动的影响。该解耦分析有助于识别各参数的主导作用。解耦分析中的参数设置见表3。

**表3. 单参数分析的参数设置**

| 变化参数 | 水平 | Asf (m?1) | K (m2) | κs (m?1) |
|---|---|---|---|---|
| 基线 | — | 1303.6 | 6.239×10?? | 498.2 |
| 比表面积 Asf | 低 | 990.7 | 6.239×10?? | 498.2 |
| | 高 | 1969.4 | 2.734×10?? | 498.2 |
| 渗透率 K | 低 | 1303.6 | 1.080×10?? | 498.2 |
| | 高 | 1303.6 | 2.734×10?? | 498.2 |
| 吸收系数 κs | 低 | 1303.6 | 6.239×10?? | 352.2 |
| | 高 | 1303.6 | 6.239×10?? | 782.6 |

图20展示了比表面积Asf、渗透率K和有效吸收系数κs独立变化时表面平均温度和蒸汽质量流率的瞬态响应。在所有情况下，表面温度在初始加热阶段快速上升，在约15–25 s达到峰值，随后随着冷却剂汽化和蒸汽排出的发展而逐渐降低。相应地，蒸汽质量流率在冷却启动后快速上升，并在准稳态阶段在约150 g/(s·m2)附近波动。与渗透率和有效吸收系数相比，比表面积对温度峰值、冷却启动和蒸汽波动的影响更为显著。

图20. 孔结构参数独立变化下表面平均物理量的瞬态响应

表4列出了孔结构参数独立变化下的冷却性能特征。其中Tmax表示透平冷却过程中的峰值表面温度；tc为冷却启动时间，定义为表面平均温度达到最大值并随后进入持续下降阶段的时间点；σm为准稳态期间蒸汽质量流率波动的标准差，使用70至100 s的数据计算得出；δl为液池厚度，由液相体积分数αl=0.5对应的中心线位置确定。比表面积对冷却性能指标的影响最为显著。当Asf从990.7增大到1969.4 m?1时，Tmax从1421.05 K升至1478.44 K，tc从16.16 s延迟至23.79 s，σm从1.463 g/(s·m2)增至4.883 g/(s·m2)，δl从2.44 mm增至2.97 mm。这些结果表明，较大的比表面积加剧了表面过热并增强了两相流的波动。相比之下，渗透率的增大通常降低峰值温度和蒸汽波动，说明较高的渗透率有利于蒸汽排出并降低两相传输阻力，尽管其对液池厚度的影响有限。κs的增大同时提高了峰值温度和液池厚度，表明更强的辐射吸收改变了多孔结构内热能量的空间沉积分布。

**表4. 孔结构参数独立变化下的特征冷却性能指标**

| 参数 | Tmax (K) | tc (s) | σm (g/(s·m2)) | δl (mm) |
|---|---|---|---|---|
| 基线 | 1432.33 | 20.30 | 3.842 | 2.79 |
| Asf=990.7 m?1 | 1421.05 | 16.16 | 1.463 | 2.44 |
| Asf=1969.4 m?1 | 1478.44 | 23.79 | 4.883 | 2.97 |
| K=1.080×10?? m2 | 1430.88 | 18.00 | 2.905 | 2.71 |
| K=2.734×10?? m2 | 1451.32 | 18.60 | 3.431 | 2.75 |
| κs=352.2 m?1 | 1427.85 | 16.74 | 2.537 | 2.56 |
| κs=782.6 m?1 | 1463.91 | 19.02 | 2.707 | 2.87 |

值得注意的是，比表面积的增大通常会扩大单位体积内的固-液接触面积并增强界面传热，按照常理预期应降低固体温度。然而，本文结果表明较大的比表面积反而导致更高的表面温度和更厚的液池。这种看似矛盾的现象表明，当辐射、相变和两相传输强烈耦合时，增强界面传热并不一定能改善整体冷却性能。

其潜在机制可能与辐射能量沉积的重新分配、液体滞留以及蒸汽排出受限有关，下文将通过固体温度、体积传热、液体体积分数和入射辐射的分布对这些机制进行进一步分析。

为了解释比表面积增大所引起的反直觉行为，图21展示了沿y方向固体温度、体积传热、液相分数和入射辐射的截面平均分布。如图21(a)(b)所示，较大的Asf导致液池区和两相区内的截面平均固体温度降低，并增强了两相区内固骨架与孔隙流体之间的热交换。由于上游区域传热更加剧烈，沿流动方向固-液温差进一步减小，从而降低了下游两相区内的液体蒸发强度。图21(c)表明，增大Asf显著提高了y = 3–5 mm范围内液相分数。较高的液相分数增强了有效辐射吸收系数，导致入射辐射降低，如图21(d)所示。此外，在较大比表面积条件下，两相区内增加的液相分数增强了局部辐射吸收，并使有效吸收系数对液相分数的瞬时变化更加敏感。因此，液相分数即使出现微小波动，也可能引起辐射能量沉积和局部相变速率的显著变化，从而改变蒸汽生成量，进一步放大相分布的波动。辐射吸收与相变之间的这种反馈最终导致蒸汽质量流量的振荡更加剧烈。

**5.3 非均匀热流密度下的特征**

在实际热防护系统中，壁面热负荷通常表现出显著的空间非均匀性。与均匀热流密度相比，非均匀热流密度不仅会导致局部热集中，还可能导致相变区域分布失衡以及冷却能力的局部不匹配。基于此，本节研究了非均匀热流密度下多孔结构内部的相变冷却，并选用孔隙率参数PPI = 40的结构B进行后续分析。非均匀热流密度的分布如图22所示。

图23(a)和(b)展示了100 s时的固体温度分布。在非均匀热流密度条件下，固体温度场沿x方向呈现出显著的温度梯度。此时，当结构受到梯形非均匀热流密度加热时，峰值固体温度达到1526.1 K，显著高于幂律热流密度分布下的值（1387.6 K），也高于均匀热流密度条件下的值（1411.0 K）。这些结果表明，即使在相同的峰值热流密度下，热流密度分布形式对冷却介质的冷却性能具有显著影响，从而导致结构峰值温度出现较大差异。冷却介质温度分布如图23(c)和(d)所示。在非均匀加热条件下，高热流侧两相区的厚度明显减小，表明该区域的冷却介质消耗更为迅速。这一现象在梯形热流密度分布下尤为明显。与幂律热流密度相比，高热流侧两相区厚度进一步减小。这表明在梯形热流密度下，冷却介质的冷却能力未能得到有效且均衡的重新分配。

图23(e)和(f)揭示了一个关键现象：在非均匀热流密度加热条件下，相变产生的蒸汽倾向于从低热流侧排出，同时将两相区内的液相冷却介质向低热流侧输送。结果是高热流侧的液态冷却介质量减少，两相区变薄，局部冷却能力减弱。这导致冷却能力与热流密度分布之间出现不匹配，尤其在高热流侧，冷却效果无法充分匹配局部热负荷。Dong和Wang [27] 也观察并讨论了类似现象。

在非均匀热流密度条件下，相变产生的蒸汽倾向于从高热流侧向低热流侧迁移并从低热流侧排出。这是因为高热流侧局部蒸发更为剧烈，导致更高的蒸汽生成速率。然而，多孔结构内的流动阻力限制了蒸汽的及时排出，促使蒸汽局部积聚，从而产生较高的局部压力。相比之下，低热流侧的蒸发较弱，局部压力也相应较低。因此，如图24所示，结构内部从高热流侧到低热流侧建立了侧向压力梯度，驱动蒸汽优先从低热流侧排出。其后果是结构内液态冷却介质的分布进一步失衡，高热流侧液体供应不足，两相区厚度减小，最终导致局部冷却性能恶化。

图25(a)展示了100 s时受热面上的冷却介质质量流量分布。在100 s时，不同热流密度分布形式下，受热面上的冷却介质质量流量呈现出显著不同的空间分布。在均匀热流密度下，质量流量沿x方向几乎保持恒定。相比之下，在非均匀热流密度下，质量流量变得强烈非均匀，表明冷却介质在多孔隙结构内发生了显著重新分配。对于梯形热流密度，质量流量分布表现出最强的空间集中性，高热流侧出现数个尖锐的局部峰值，向低热流侧逐渐下降，最终达到65 g/(s·m2)的值。这表明冷却介质的蒸发和消耗表现出强烈的局部特性，导致冷却能力与热负荷之间出现不匹配。幂律热流密度也产生了非均匀的质量流量分布，但变化相对平缓。此外，如图25(b)所示的冷却介质瞬时质量流量表明，与均匀热流密度相比，非均匀加热条件下蒸汽质量流量的波动更大。

总体而言，质量流量分布进一步证实了不均匀加热不仅改变了温度场和两相分布，还显著改变了受热面上的蒸汽质量流量分布，导致冷却介质流量与冷却需求之间出现不匹配。

**6. 结论**

本文建立了物理模型，并在OpenFOAM中开发了数值计算代码，研究了透气冷却应用中多孔结构内液态冷却介质相变冷却的瞬态传热传质特性。在此基础上，系统研究了热辐射、孔隙参数PPI以及非均匀热流密度分布对多孔结构内相变冷却过程的影响，得出以下有意义的结论。

辐射传热是高温透气冷却过程中的重要机制。它对局部温升、液池区域位置、冷却启动过程和峰值温度有显著影响。随着吸收系数的增加，辐射穿透深度减小，辐射能量更容易集中在结构B的上部区域。这显著影响了固体温度分布和两相区的演化特性，导致冷却启动延迟、峰值温度升高以及结构内部温度梯度增大。

在独立变化的孔隙相关参数中，比表面积对瞬态冷却行为的影响最为显著。将Asf从990.7增大到1969.4 m?1，峰值温度从1421.05 K升高到1478.44 K，冷却启动时间从16.16 s延迟到23.79 s，液池厚度从2.44 mm增加到2.97 mm，蒸汽质量流量波动的标准差从1.463 g/(s·m2)放大到4.883 g/(s·m2)。虽然较大的比表面积增强了固-流体传热，但也增加了两相区内的液体滞留，并降低了下游固-液温差，抑制了蒸发，增加了两相区内的液相分数。较高的液相分数增强了辐射吸收并减弱了下游辐射传输。其波动进一步改变了辐射能量沉积和局部蒸发，从而加剧了蒸汽生成和排出振荡。总体而言，冷却行为强烈受固-流体传热、辐射传输和相变之间耦合关系的控制。

在非均匀热流密度条件下，高热流侧相比低热流侧局部蒸发更为剧烈，导致更高的蒸汽生成速率。然而，由于多孔结构内的流动阻力，蒸汽无法及时排出并倾向于局部积聚，从而产生较高的局部压力。因此，结构内部从高热流侧到低热流侧建立了侧向压力梯度。该压力梯度驱动蒸汽向低热流侧迁移并排出。这导致了冷却能力与施加热负荷之间的显著不匹配。

**作者贡献声明**

Xin Lin Xia：写作–审阅与编辑、调查研究、资金获取、概念构思。

Wei Du：写作–审阅与编辑、写作–原稿撰写、软件、调查研究。

Qin Yang Wang：资源、项目管理。

Xian Ying Chu：资源、项目管理。

Xiao Lei Li：资金获取。

Wei Du | Qin Yang Wang | Xin Lin Xia | Xian Ying Chu | Xiao Lei Li

哈尔滨工业大学能源科学与工程学院，工信部航空航天热物理与涡轮动力重点实验室，中国黑龙江省哈尔滨市 150001