吴 波，张丰华，田 沣，醋强一

（中国航空工业集团有限公司西安航空计算技术研究所，陕西西安710068）

航空电子模块是战斗机、空空导弹等航空飞行器中实现信息测量、采集、传输、处理、监控及显示功能的基本载体。新一代战斗机、空空导弹等的性能更加强大［1-2］，其航空电子模块的功能越来越复杂，总功耗也越来越大。而“小、低、轻”是对航空电子模块的必然要求，这使得航空电子模块的热管理问题越来越突出：一方面飞行器的性能增强导致功耗增大；另一方面其自身的空间尺寸和质量受到严格的限制。因此，通过引入外部工质实施对流冷却的方案在很多情况下并不可行，例如远端安装的航空电子模块不便于流体管路的布置。相变储热器是一类依靠相变材料吸热熔化来储存热量的装置，可与航空电子模块集成在一起，以便实施热管理［3］。当航空电子模块具有短时工作特点时，相变储热器更能发挥作用。

工程师在设计储热器时需要关注的核心问题是：根据热源芯片的许用温度和功耗来判断储热器从航空电子模块开始加电运行到热源芯片达到许用温度的时间，即储热器的工作时间。采用数值传热学方法及CFD（computational fluid dynamics，计算流体动力学）软件可以对储热器的工作时间进行数值仿真计算。但是，相变传热过程的精确数值模拟较为复杂，不适合应用于需要迅速响应的工程设计中，也不便在非传热学背景的结构设计工程师中推广。因此，为了满足工程设计中对设计工具的“易用”“速用”需求，笔者拟构建储热器工作时间与功耗、设定温度的数学关系式。

目前，已有大量关于类似相变储热器结构的基础研究和应用研究［4-10］，其可为本文的研究提供重要参考。但是，这些研究未综合考虑航空电子模块储热器实际工作时与自然对流和热辐射的耦合作用，也缺乏在航空电子模块热管理工程背景下对储热器工作时间的预测。

基于此，笔者拟开展以下工作。首先，综合考虑航空电子模块储热器在实际工作时其顶面的自然对流和热辐射，运用传热学理论构建储热器热模型。然后，运用焓-多孔介质法对储热器进行数值模拟，并以顶面绝热条件下的储热器为参照，分析储热器在工作过程中的温升规律及其相变材料的熔化进程。接着，根据已有的相变储热器研究结果，结合数值模拟结果，分析储热器工作时间与功耗、设定温度之间的依变关系。最后，运用二元线性回归方法拟合得到储热器工作时间的预测公式，并通过独立的数值模拟数据集来验证该预测公式的准确性。

1 航空电子模块储热器热模型

航空电子模块储热器是通过焊接密封的方式将固-液相变材料封装在金属腔体内部的。金属腔体内部有1组等间距的肋片，用于增强腔体内部的换热。如图1所示，航空电子模块由印制电路板（printed circuit board，PCB）和电子芯片组成，并通过紧固件与储热器连接。电子芯片工作时产生的热量通过热传导的方式传输给储热器，储热器内部的相变材料吸收热量。当相变材料达到熔点时，其温度几乎保持不变，从而实现对电子芯片的温度控制。储热器内部每2个肋片对称面之间的部分为储热器的最小单元，即图1 中虚线框内部分。假设电子芯片的热量均匀地分布在储热器底面上，则该最小单元可视作表征储热器储热过程的特征单元体，通过研究该特征单元体工作时间的特性，即可得到整个储热器工作时间的特性。储热器单元体的几何模型如图1下方所示，其几何尺寸如表1所示，其中L为储热器的总长度，a和b分别为储热器单元体左、右侧和上、下侧金属壁的厚度，H和l 分别为储热器单元体金属腔体的高度和宽度。储热器单元体金属壁面的材料为铝6061，内部相变材料为RT80，其热物性参数详见文献［11］。

图1 航空电子模块储热器热模型Fig.1 Thermal model of heat storage device in avionics module

表1 储热器几何尺寸Table 1 Geometric dimensions of heat storage device

储热器单元体的底面接收由电子芯片产生的均匀、热流密度恒定的热量，顶面直接暴露在外部环境中，通过自然对流和热辐射向环境散热，两者的热量差以相变材料潜热的形式储存在储热器中。假设储热器的总长度和总宽度均为L，则其底面的输入热流密度qin可表示为：

式中：Q 为航空电子模块的总功耗，本文考察4 种典型功耗，即Q=12，15，21和30 W。

储热器顶面的输出热流密度qout由自然对流热流密度qnc和热辐射热流密度qrad两部分组成，即：

其中：

式中：h为储热器顶面的自然对流换热系数；Ttop为储热器顶面的平均温度；Tamb为环境温度，本文取Tamb=333.15 K；λa为环境温度下空气的导热系数；Nu、Gr和Pr分别为努赛尔数、格拉晓夫数和普朗特数，其具体计算公式详见文献［12］，计算这几个无量纲准则数时采用环境温度下空气的物性参数，储热器特征尺寸均取L。

根据斯蒂芬-玻尔兹曼定律，热辐射热流密度qrad可表示为：式中：ε为储热器顶面的发射率，本文取ε=0.8。

储热器稳态工作时，其底面的平均温度与电子芯片的温度直接相关。在实际应用中，由于热阻的作用，前者略低于后者。考虑到安全余量，可直接将储热器底面的平均温度视作温度控制的目标温度。本文定义从电子芯片开始产热到储热器底面平均温度达到设定温度的时间为储热器的工作时间。

2 储热器工作时间预测的研究方法

2.1 数值模拟方法

本文采用焓-多孔介质法来模拟储热器单元体内部相变材料的二维瞬态熔化过程。在模拟相变材料的熔化过程时，考虑液相相变材料的自然对流流动，并假定该流动为层流。假设储热器内部相变材料的固相与液相之间不存在密度差异，其浮升力项中的密度按照包辛尼斯假设处理；金属壁面内导热过程为二维瞬态纯导热过程；忽略金属壁面与相变材料之间的热阻。储热器中相变材料区域的控制方程如表2所示。

表2 储热器中相变材料区域的控制方程Table 2 Governing equations of phase change material zone in heat storage device

表2 中：u、v 分别为x、y 方向的速度；p 为压力；T为温度；t为时间；ρP、μP、λP和CP分别为相变材料的密度、黏度、热导率和比热容；Su、Sv和ST分别为各控制方程的源项，其具体表达式如下：

式中：g为重力加速度；β、E分别为相变材料的热膨胀系数和焓；A为多孔系数；Tref为相变材料的参考温度，本文中Tref的值为相变材料的初始温度；ρPref为参考温度下相变材料的密度［13］。

设置储热器单元体的边界条件为：左、右两侧金属壁面为对称边界；底面为恒热流密度边界，热流密度大小由式（1）求得；当考虑自然对流和热辐射时，顶面为变热流密度边界，热流密度由式（2）至式（5）求得；当不考虑自然对流和热辐射时，顶面为绝热边界。在数值模拟求解过程中，设相变材料的初始温度Ti=343.15 K。

基于上述控制方程和边界条件，在Fluent计算平台上对储热器单元体内部相变材料的二维瞬态熔化过程进行数值模拟。在进行数值模拟前，先对采用的数值模拟方法的准确性进行验证，验证方法为文献［14］中的方法，即用具有温度场理论解的斯蒂芬问题来验证本文采用的焓-多孔介质法的准确性。利用焓-多孔介质法求解经典的斯蒂芬问题，计算域大小为0.01 m×0.10 m。对同一时刻计算域长度方向温度的数值模拟值与理论值［15］进行对比，结果如图2 所示。从图2 中可以看出，数值模拟值与理论值较吻合，两者的最大误差为0.54%，说明本文采用的数值模拟方法具有一定的准确性。

图2 焓-多孔介质法的准确性验证结果Fig.2 Ⅴerification results of accuracy of enthalpy-porous media method

2.2 线性回归方法

航空电子模块储热器的工作时间会随功耗和设定温度的变化而变化，可利用线性回归方法对其进行预测。

当因变量Y随2个自变量x1、x2变化时，二元线性回归预测方程可表示为：

式中：Yi、x1i和x2i分别为Y、x1和x2的观测数据。

为了考察二元线性回归方法的预测精度，采用拟合优度R2进行评价，其表达式如下：

3 结果分析

3.1 储热器温升及其相变材料熔化进程的规律

在不考虑储热器顶面自然对流和热辐射（即顶面绝热）时，4 种功耗下储热器底面平均温度随时间变化的曲线如图3所示。从图3中可以看出，随着功耗的减小，相变材料熔化过程持续的时间延长，且在其熔化过程中储热器底面平均温度的变化越平缓。在不同功耗下，在相变材料开始熔化前和熔化结束后，储热器底面的平均温度均线性升高。

图3 顶面绝热时不同功耗下储热器底面的温升曲线Fig.3 Temperature rise curves of bottom surface of heat storage device under different power consumptions with top surface insulation

在考虑储热器顶面自然对流和热辐射（即顶面散热）时，4 种功耗下储热器底面平均温度随时间变化的曲线如图4所示。对比图3和图4可知，在相同功耗下，考虑顶面散热时，储热器底面平均温度的上升更为缓慢；在相变材料熔化结束后，储热器底面平均温度的升高偏离线性变化，且功耗越小偏离程度越大。这是因为功耗越小，储热器顶面通过自然对流和热辐射所散的热量的占比越高，且自然对流和热辐射并不恒定，随着温度的升高，通过自然对流和热辐射所散的热量越大，从而导致储热器底面的平均温度上升越缓慢。

图4 顶面散热时不同功耗下储热器底面的温升曲线Fig.4 Temperature rise curves of bottom surface of heat storage device under different power consumptions with top surface heat dissipation

为了深入了解相变材料的熔化进程，进一步考察2种情况（不考虑/考虑储热器顶面自然对流和热辐射）下储热器内部相变材料的熔化过程。提取不同功耗下相变材料熔化过程中液相组分含量F的瞬态变化曲线，即相变材料熔化进程曲线，结果如图5所示。其中液相组分含量F为相变材料熔化过程中液相体积与总体积的比值，F=0表示相变材料未开始熔化，F=1表示相变材料熔化结束。从图5中可以看出，考虑顶面自然对流和热辐射时，储热器内部相变材料熔化进程变慢，但当功耗为30 W时，储热器内部相变材料熔化进程的变慢程度并不明显；当功耗减小为12 W 时，储热器内部相变材料熔化进程的变慢程度越来越显著。

为了得到储热器工作时间的普遍性规律，本文对时间、功耗和温度进行无量纲化。采用傅里叶数Fo作为时间的无量纲量，为：

式中：α为相变材料的热扩散系数。

采用相变材料的初始温度Ti及熔点温度Tm（Tm=353.15 K）对温度T进行无量纲化。定义无量纲温度θ为：

图5 不同功耗下相变材料的熔化进程曲线Fig.5 Melting process curves of phase change material under different power consumptions

根据图4所示的曲线，获取考虑顶面自然对流和热辐射时储热器在不同功耗及设定温度下的工作时间。为了给下文的预测提供数据，从图4的曲线中采样并构建2组数据：一组包含20个数据，作为建模集，采样功耗分别为30，21，15和12 W，采样温度分别为365，370，375，380和385 K；另一组包含16个数据，作为验证集，采样功耗分别为30，21，15 和12 W，采样温度分别为367.5，372.5，377.5和382.5 K。

3.2 储热器工作时间与功耗、设定温度的关系

从上文结果可知，储热器的工作时间除了与功耗有关，还与设定温度有关。为了准确预测储热器的工作时间，需分析其与功耗、设定温度的关系。

3.2.1 工作时间与功耗的关系

从文献［16-18］中可知，在不考虑储热器顶面自然对流和热辐射时，在不同功耗下储热器内部相变材料熔化过程中，采用斯蒂芬数Ste 和瑞利数Ra 修正后，无量纲工作时间Fo与液相组分含量F 之间在不同功耗下具有相同的函数关系，即：

从文献［16］中可知，斯蒂芬数Ste和瑞利数Ra与无量纲功耗ω 均呈正比，而α为常数。由此可见，当储热器内部相变材料的液相组分含量F确定时，由于F 与无量纲功耗ω 无关，则无量纲工作时间Fo 与无量纲功耗ω的（1+α）次方呈反比。换言之，在不考虑顶面自然对流和热辐射时，储热器的工作时间与功耗呈指数关系。

根据上述分析可知，在考虑顶面自然对流和热辐射的情况下，基于建模集中的20个数据，绘制无量纲工作时间Fo与无量纲功耗ω在单对数坐标（纵坐标为以10为底的对数坐标，横坐标为线性坐标）下的依变关系曲线，结果如图6所示。

图6 单对数坐标下Fo随ω的变化规律Fig.6 Ⅴariation law of Fo with ω in single logarithmic coordinate

由图6可知，在单对数坐标下，无量纲工作时间Fo与无量纲功耗ω之间不呈线性关系。这是因为在不考虑顶面自然对流和热辐射时，输入储热器单元体的净热量为定值；而在考虑顶面自然对流和热辐射时，输入储热器单元体的净热量随温度的升高而加速减小。因此，在图6的基础上，对无量纲功耗ω 也取对数，结果如图7所示。

图7 双对数坐标下Fo随ω的变化规律Fig.7 Ⅴariation law of Fo with ω in double logarithmic coordinate

从图7中可以发现，在双对数坐标下，无量纲工作时间Fo与无量纲功耗ω呈线性关系，且无量纲设定温度θ 取不同值时，4 条变化曲线几乎彼此平行。图7所示的无量纲工作时间Fo与无量纲功耗ω之间的线性关系可近似用微分方程表示为：

3.2.2 工作时间与设定温度的关系

从图4中可看出，储热器内部相变材料熔化结束后，储热器底面平均温度随时间的变化曲线近似为对数曲线。因此，在考虑顶面自然对流和热辐射的情况下，基于建模集中的20个数据，绘制无量纲工作时间Fo与无量纲设定温度θ在单对数坐标下的依变关系曲线，结果如图8所示。

图8 单对数坐标下Fo随θ的变化规律Fig.8 Ⅴariation law of Fo with θ in single logarithmic coordinate

从图8中可看出，在单对数坐标下，无量纲工作时间Fo 与无量纲设定温度θ 几乎呈线性关系，且无量纲功耗ω取不同值时，各条曲线几乎平行，则无量纲工作时间Fo与无量纲设定温度θ之间的线性关系可近似用微分方程表示为：

3.3 储热器工作时间的回归预测

由此可见，可利用二元线性回归方法对储热器工作时间进行预测。基于建模集中的20个数据，求解方程组（10），得到b0、b1和b2。本文运用R语言中的二元线性回归函数来求解b0、b1和b2，最终得到的拟合关系式为：

通过二元线性回归拟合，得到建模集中20 个数据对应的lg Fo 的拟合值，并与对应的数值模拟值进行对比，结果如图9 所示。从图9中可以看出，20个lg Fo的拟合值非常接近其数值模拟值，且最大误差在±4%以内，拟合优度R2=0.993。

图9 建模集中lg Fo的拟合值与数值模拟值对比Fig.9 Comparison of fitted values and numerical simulation values of lg Fo in the discovery set

为了验证式（23）的预测精度，利用由16 个数据组成的验证集进行验证。利用式（23）计算验证集中16个数据对应的lg Fo预测值，并与数值模拟值进行对比，结果如图10所示。从图10中可看出，式（23）的预测精度非常高，最大误差为±3%。

图10 验证集中lg Fo的预测值与数值模拟值对比Fig.10 Comparison of predicted values and numerical simulation values of lg Fo in the validation set

4 结束语

基于建立的航空电子模块储热器在实际工作环境下的热模型，利用数值模拟方法分析发现，储热器工作时间与功耗、设定温度之间在对数坐标下呈线性关系，进而利用二元线性回归法拟合得到储热器工作时间的预测公式。

基于拟合得到的储热器工作时间预测公式，设计人员可方便地计算出某一功耗和设定温度下航空电子模块储热器的工作时间，这对于弹载类短时工作的航空电子模块的热管理设计具有重要意义，可避免工程设计中不必要的、复杂的、非线性且多散热路径耦合的相变传热数值模拟。该预测公式可作为航空电子模块储热器快速评估和敏捷开发的工具。但是，需要注意的是，该预测公式是在给定储热器尺寸下得到的，并未考虑尺寸变化的影响，后续研究中可对此进行完善。