深冷环路热管瞬态数值仿真研究

何发龙，杜王芳，赵建福，3*，陈跃勇，何 江，张红星，苗建印

(1.中国科学院力学研究所微重力重点实验室,北京100190；2.空间热控技术北京市重点实验室,北京100094；3.中国科学院大学工程科学学院,北京100049；4.中国空间技术研究院北京空间飞行器总体设计部,北京100094)

1 引言

深冷环路热管(Cryogenic Loop Heat Pipe,CLHP)作为空间红外探测系统关键组件,可以有效隔离制冷器机械振动及电磁干扰对光学探测仪器成像造成的影响［1］,在工程实践中有着重要应用。目前,以氦作为工质,可完成CLHP超临界启动,并实现3～5 K系统运行温度条件下10～140 mW热量传输［2］。此外,对于其他深冷工质的CLHP也有相应实验研究,均可实现系统超临界启动并维持一定温度范围内稳定运行,如氮(80～110 K)［3］、氖(30～40 K)［4-5］、 氢 (20～30 K)［6］。 目 前, 对 于CLHP的实验研究主要围绕系统的超临界启动和稳定运行性能分析。

然而,CLHP技术研发和实际运控面临的最大技术挑战在于缺乏对CLHP内部多相热流体动力学行为特性及其机理的可靠认知,以及对相关现象的可靠、合理的数学描述。尤其是CLHP运行过程中各部组件内部流动特征(工质分布特征、单相流态和/或气液两相流型、能量平衡与压力平衡特性等),以及CLHP内部流体工质分布、热力状态变化及其流动与传热规律等。

相比实验研究,理论分析和数值模拟则是定性揭示或解释CLHP传热机制的重要手段。CLHP作为一种特殊结构形式的环路热管,用于模拟分析的数值模型与环路热管相近,可分为稳态模型和瞬态模型。Kaya［7］基于环路热管各部件的能量平衡建立了一维稳态分析模型,并模拟分析出环路热管的变热导区性能。柏立战等［8］考虑到冷凝器内两相流型的因素,基于环状流型,建立了一种稳态分析模型,可用于确定恒定热导运行模式下储液器内液相体积增加造成的冷凝器内部两相区的长度缩减。随后将该模型扩展至CLHP的稳态分析［9］,并进行了氮工质的 CLHP稳态运行参数分析。稳态分析模型可用于分析CLHP系统的稳定运行性能和参数影响分析,但无法表征系统内部各组件的流动特征。因而,Kaya［10］采用交错格点法建立一种环路热管瞬态模型,可用于预测系统瞬态响应,但该模型中未考虑到对蒸气槽道和蒸发器内核的工质状态变化,因而毛细芯内部的流动换热参数难以确定,特别是在低热负荷条件下,蒸发器逆向漏热量需要通过实验值拟合确定。为了增强模型适用性就需要减少模型中的实验拟合参数,Nishikawara［11］在该模型的基础上,将蒸气槽道和蒸发器内核作为节点建立数值模型,但模型中的动量方程是以稳态形式给出,因此是一种近似稳态模型。柏立战［12］建立一种CLHP超临界启动模型,并用以分析了工质充装压力、寄生漏热等因素对于超临界启动的影响,模型中为避免超临界启动过程中工质温度大幅变化造成的热物性变化,从而采用温度拟合的多项式来计算工质热物性,虽然简化了模型计算,但与工质实际的热物性状态差距较大,从而降低了模型的计算精度。此外,模型并不能完整实现CLHP超临界启动的全过程模拟,需要根据不同运行阶段重新设定初始条件以完成模拟。

本文基于具有辅回路的CLHP结构特征及其运行过程特点,构建一个一维CLHP瞬态数值模型,通过数值仿真分析CLHP稳态运行中流体工质的分布及其热力学状态变化路径,并据此对CLHP稳定运行所必须满足的条件进行分析。

2 瞬态数值模型

2.1 基本假设

本文研究对象是图1所示的具有辅回路的CLHP系统［13］,其瞬态仿真模型的构建采用节点网络法,并基于组件内部工质分布、热力状态变化及流动传热特性［14］,采用了如下基本假设：

1)工质节点处流体总是处于局部热力学平衡状态,且满足吉布斯相律,即单组分工质只有2个独立的热力学状态参数(温度和密度),其它热力学状态参数采用查表计算；

2)毛细芯骨架温度和毛细芯内部工质温度相同,忽略气-固和固-液界面接触热阻以及固体骨架热惯性的影响。

2.2 控制方程

CLHP系统瞬态仿真模型包括相互耦合的热力学模型和动力学模型,其中,工质节点的热力学模型采用控制容积方法,并基于质量守恒定律和能量守恒定律建立；而沿流动路径(即相邻工质节点间的管段)工质的动力学模型同样采用控制容积方法,并基于动量守恒定律建立。至于壁面节点只需要基于能量守恒定律建立描述其温度变化的热力学模型即可。

工质节点i质量方程如式(1)所示：

其中,mi为工质节点i的质量,m·ik为邻近工质节点k与工质节点i界面处的质量流率。

工质节点i能量方程如式(2)所示：

其中,hi为工质节点i的焓值,m·ikhik为邻近工质节点k与工质节点i界面处的焓流率(这里,hik的取值依赖于m·ik的方向,即总取流动上游的工质节点比焓),Qim为工质节点i与相邻壁面节点m之间的换热量,Wi为工质节点i因压力变化引起的压缩功。

流动路径ik(即相邻工质节点i和k之间流动管段)上工质的动量方程如式(3)所示：

其中,Pup和Pdown分别表示该路径上游和下游工质节点的压力,ΔPa、ΔPf和ΔPg分别表示该路径上的加速压降、摩擦压降和重位压降,Leff,ik和Aeff,ik分别表示该路径管段的等效长度和等效流通面积。

壁面节点j能量方程如式(4)所示：

其中,mj和Cj分别表示壁面节点j的质量和比热,Qjn为该壁面节点与相邻的其它壁面节点、工质节点及外部环境或施加的热源/热沉之间的换热量。

在蒸发器工作时,需引入气液相变界面节点,其控制方程如式(5)所示：

工质与壁面的流动换热系数为式(6)：

其中,h为对流换热系数,kf为工质的导热系数,Nu为努赛尔数。流动处于层流(Re≤1960)时,Nu＝3.66；否则,流动为湍流(Re＞6420),Nu可由Dittus-Boelter公式计算,如式(7)所示：

其中,工质被冷却时n＝0.3,工质被加热时n＝0.4。而工质流动路径上的压力损失如式(8)所示：

其中,摩擦阻力系数 f由 Churchill［15］公式确定。

毛细芯内的流动阻力为式(9)：

其中,Kp为毛细芯渗透率,m为工质动力粘度,Ro,wi为毛细芯外径,Ri,wi为毛细芯内径,Lwi为毛细芯长度。

上述微分-代数方程构成CLHP瞬态数值模型。结合流体工质物性参数库,以及适当的边界条件(一般按照实际热边界状态予以设置)和初始条件(如超临界启动计算时可设置为均匀的超临界状态),即构成对CLHP系统的完整描述。

2.3 模型求解

对上述CLHP瞬态仿真模型,本文采用隐式变步长微分代数计算法DASSL［16］进行求解,方程中的导数采用1～5阶向后差分公式计算,每个求解步长内,求解器根据设定误差限自动选择合适的阶数和步长,在维持计算稳定性的同时,减少迭代次数。

2.4 仿真参数设定

为了与测试结果进行比较,本文仿真采用了与何江等［4］、郭元东等［5］试验系统相同的参数：工质氖,工作温区30～40 K；除毛细芯体为纯镍粉烧结外,回路中其他部件均采用不锈钢。结构参数见表1。

表1 CLHP的基本参数Table 1 Basic parameters of CLHP

本文计算中,将储液器、蒸发器、蒸气腔、加热源、储气室作为单一节点,主辅回路传输管线及冷凝器划分19个节点,共计29个节点,模拟工况设置如下：

系统充装压力 Psys＝2.1 MPa,环境温度Tam＝290 K,热沉温度Tsink＝35 K,热沉与冷凝器及次储液器之间的热导hc,s＝240 W/K,暂不考虑环境漏热影响,同时,为简化起见也未考虑重力的影响。初始条件根据系统可能的物理状态人为设定。需要说明的是,尽管所设定的初始状态并不能保证在流动与传热细节上完全符合真实物理约束,但其中非物理真实引起的差异会在求解过程中自动消除,对计算结果没有影响。

郭元东等［5］试验研究了辅蒸发器热载荷变化对系统传输能力的影响,发现增大辅蒸发器热载荷会减弱系统热传输能力,当辅蒸发器热载荷大于1.5 W时,这种恶化趋势更明显。另一方面,辅蒸发器热载荷低于1.5 W时,系统超临界启动难以成功。因而,本文计算中辅蒸发器热载荷取值为QSE＝1.5 W。主蒸发器热载荷QPE采用阶梯上升方式给定,分别为0.5 W、1 W、1.5 W、

2 W、3 W和3.5 W,以确定其传热性能曲线。

3 仿真结果与讨论

3.1 蒸发器温度响应

图2和图3给出了数值仿真得到的主蒸发器壁面温度TPE,W、辅蒸发器壁面温度TSE,W和主蒸发器槽道内工质温度TPE,V、辅蒸发器槽道内工质温度TSE,V及储气室压力Psys随时间的变化情况。仿真计算中,由于人为设置的初始条件难以避免非物理性偏差的存在,仿真过程需要额外时间消除其影响,这导致计算开始时系统趋于稳定需要更多时间,相比于后续热载荷提升时要明显缓慢得多。

在本文研究参数范围内,主蒸发器热载荷增加导致主蒸发器壁面温度依次升高,系统压力仅在1.5 W时略有下降,这主要是系统内气液质量比减小造成的,随着热负荷的增加,产生的蒸气量增加,在管径不变的条件下,冷凝器中的工质流动状态向湍流过度增加了冷凝器的冷凝效率,更多的蒸气冷凝为液相增加了系统中的液相比重。但压力总体上呈现依次提高的趋势,而基于局部平衡假设,工质温度与系统压力变化同步,趋势相同。由于辅蒸发器热载荷维持恒定不变,辅蒸发器壁面温度和工质温度之差维持恒定,但同样存在随工质温度变化而依次提高的趋势。

3.2 系统传热性能

图4给出了基于主蒸发器壁面温度TPE,W和主蒸发器热载荷QPE的CLHP系统传热性能曲线。不考虑环境漏热及重力的影响,CLHP工作温度与主蒸发器热载荷近似呈线性关系,只在低热载荷时存在一定的偏离,工作温度会略高于线性趋势预言的数值,这与郭元东［5］试验观测是一致的,后者观测到CLHP在低负荷区稳定运行时工作温度会略高于线性预测值。

3.3 变热载荷条件下的系统压降

图5 给出了系统主、辅回路总压降ΔPPL,T和ΔPSL,T随主蒸发器热载荷QPE的变化情况。当主蒸发器热载荷较小时,主回路压降低于辅回路压降；随着主蒸发器热载荷的增加,系统主、辅回路的压降也随之增加,且主回路压降增加得更快；工作温度的增加则导致表面张力的减小,因此最大毛细抽力ΔPcap,max略有下降。原则上,CLHP稳态运行时最大毛细抽力不小于任一回路压降。图5中,QPE＝3.5 W时最大毛细抽力略小于主回路压降,显然,毛细抽力的不足,会使得该状态难以稳定维持。计算中,人为保持该工况,是为了确定CLHP传热极限,即在本文仿真分析的参数条件下,CLHP传热能力将略低于3.5 W。郭元东等［5］研究了充装压力对CLHP传热能力的影响,相同条件下的数值约为3.0 W,其间差异0.5 W与系统寄生漏热量基本一致,从另一个方面验证了本文数值仿真结果的可靠性。

3.4 工质分布状态与系统热力学状态变化路径

图6 给出了主蒸发器热载荷3 W时系统稳态运行中内部工质的热力学状态变化(不包括储气室及其连接管线)。图中,点A和A*、B和B*分别对应主、辅蒸发器毛细芯内气液界面的气相侧和液相侧,其间的压力差对应着相应回路压降,稳态运行时其数值不能大于最大毛细抽吸力。由图可知,蒸发器受热导致液体工质蒸发产生饱和蒸气,蒸气在槽道和气相管线流动中继续吸收热量和/或压力下降而转变为(稳态)过热蒸气。过热蒸气进入主冷凝器内放热凝结,转变为饱和状态,经主液相管线进入主储液器及辅回路,而主、辅储液器内的饱和液体在毛细抽吸力作用下再次进入蒸发器毛细芯内,随压力的下降和毛细芯导热转化为亚稳态的过热液体,并在气液界面处再次蒸发,形成工质的主、辅回路循环。除储气室及其连接管路外,系统内部工质温度在稳态运行条件下最大变化不超过0.2 K,因此,整个系统将具有极为良好的均温特性。

图6中,含液管线内部工质总是处于饱和状态,其温度随压力的沿程下降而随之降低,其路径与工质气液饱和线相同,这明显是气液混合状态的表现。查验系统内工质分布状态(图7)可以发现,主、辅回路中冷凝器、储液器和液体管线内气相干度均处于0～1之间,确实一直处于气液两相混合状态。需要注意的是,尽管初始条件假设了干度为0的纯饱和液相状态,仿真计算中会迅速演变为气液两相混合状态。热载荷变化时,干度大小会有所调整,但总体状态并未发生变化。

4 关于CLHP使役条件的进一步讨论

Maydanik［17］系统阐述了环路热管的使役条件,即能够稳定工作的必要条件。赵建福等［14］将其推广到具有辅回路的CLHP系统,并给出了CLHP稳定工作的条件：

1)主、辅回路总压降均应不大于相应的主、辅蒸发器毛细芯提供的最大毛细压力,如式(10)所示：

式中,I＝1和2,分别代表主回路和辅回路。这与一般热管要求相同,对应于毛细极限情形。

2)主、辅蒸发器毛细芯蒸发面和相应的储液器内必须保持足够的压差和温差,且满足式(11)所示的关系：

式中,J＝a或s,分别代表主、辅蒸发器蒸气槽道饱和蒸气温度；K＝k或q,分别代表主、辅储液器中饱和蒸气温度。该条件要求确保启动时在储液器里能有足够液体工质。

3)主、辅液相回路不能出现过热状态,即式(12)：

式中,L＝h或n分别代表主、辅冷凝器内的蒸气完全凝结位置,M＝j或p则分别代表主、辅液相管路出口。该条件的基本要求是液体工质不会因压力下降和环境漏热而发生沸腾。

不过,关于使役条件(3)的一般解读,往往沿袭Maydanik［17］的说法,即主、辅液相回路必须保证足够的过冷度。但是,由图6和图7可以明确看到,“足够的过冷度”要求是对式(12)的过分解读,事实上,条件(11)、(12)并不排斥等号的存在,即饱和状态也自然是满足使役条件的。因此,图8示意出的主、辅液相回路中的过冷状态,应该并不排斥图6中的饱和变化路径。

周顺涛［18］实测了氮工质 CLHP稳态运行时系统温度分布,在较高充装压力时,主冷凝器出口温度显著低于主液相管线出口温度,后者又明显低于主储液器温度,表明在主液相管线内流动着的是干度为0的纯液相工质,且确实处于过冷状态,寄生漏热导致工质和管壁温度沿程不断升高；然而,在较低充装压力时,主冷凝器出口温度要略高于主储液器温度,呈现出和本文仿真结果相同的趋势。后者只能由气液两相混合状态来解释：主冷凝器出口压力要高于下游的主储液器,相应的饱和温度也会略高。

周顺涛［18］还发现在较高充装压力时,主液相管线出口温度存在明显毛刺现象,这被认为是因为过冷沸腾的发生引起的,即过冷工质因寄生漏热与热载荷变动等相互作用产生沸腾,管线内气泡的产生引起系统显著波动。这会使得系统在明显低于其传热极限时失稳,导致运行失败。与此相反,在较低充装压力时,主液相管线温度变化平稳,扰动被抑制。在郭元东等［5］实验中也同样观测到了类似的充装压力对氖工质CLHP稳定性的影响现象。

上述现象可以从以下事实得以解释：主冷凝器出口过冷的纯液体工质流动过程中受寄生漏热影响,会在局部达到过冷沸腾条件,鉴于CLHP内部往往充装的是高纯度深低温工质,对管壁润湿性优异,相应的沸腾起始过热度会较高,过冷沸腾发生时伴随气泡生成释放出来的扰动强度大,对系统稳定性的扰动强,进而弱化了系统稳定性。而当含液管线处于气液两相混合状态时,天然存在的气液界面,使得液相受热蒸发能够平稳发生,从而大大抑制了寄生漏热对系统稳定性的扰动。因此,若仅考虑CLHP稳定性,含液管线内工质并不宜处于过冷状态。当然,较低充装压力时,传热极限易受毛细芯供液量不足的限制,需要综合考虑各方面特性需要,开展优化设计。

5 结论

本文基于对CLHP工作原理与内部工质的热动力学规律的认识,采用节点网络法构建了具有辅回路的CLHP瞬态数值仿真模型,对稳定热边界条件下氖工质CLHP运行性能进行了非稳态仿真计算,得到了CLHP稳态工作性能曲线和系统内部工质分布状态、热力学状态变化路径及系统压降等,数值结果与试验观测吻合较好,表明所构建的CLHP瞬态数值仿真模型能够描述具有辅回路的CLHP内部流动与传热特性,可以用于系统性能的仿真分析。

基于数值仿真结果,特别是对含液管线内部工质分布状态的预测结果,本文对CLHP使役条件进行了深入分析,指出传统观点要求的“液相回路必须保证足够的过冷度”属于过度解读,其结果在CLHP无法完全隔绝寄生漏热的条件下,反而会因液体工质的过冷沸腾消弱系统稳定性。