液氮重力自循环系统的仿真及设计
2018-05-08王安良
高 天，王安良
(北京航空航天大学宇航学院， 北京 100191)
simulation
液氮重力自循环系统作为一种低温热沉系统，其工作原理如图1所示，其主要物理机制是热虹吸流和自然对流沸腾。流体在传热的作用下流动，热沉管内液氮因壁面受热而产生自然对流或部分汽化，致使上升段与下降段产生重力压差，且流动同时受摩擦阻力的限制，流动达到稳定时重力压差与系统总阻力达到动态平衡。
国内，重力自循环系统直接应用于工程领域的研究，工作较少且集中在结构设计与实验方面[1]。国外，Podkorytov等[2]实验研究了毛细力对液氮自然对流循环过程中传热的影响。韩国Choi等[3]、Kim和Chang[4]以及Chang[5]课题组把液氮重力自循环制冷方法用于高温超导研究领域，他们的研究结果表明，液氮重力自循环系统可用于高温超导(high temperature superconductor,HTS)领域，如果液氮有一定的过冷度，效果更好；随着表面热流密度的升高，被冷却壁面温度也在升高。日本中部大学联合美国MIT等多家单位[6-9]针对过冷液氮重力自循环系统在HTS领域的应用开展了实验与理论研究工作。
法国的Baudouy[10]开展了低热流密度下液氮自然循环系统的实验研究，并采用Gr数，Nu数分析了两相流区段的换热特性，并给出了一个修正的计算自然对流关联式的公式
(1)
式中，NuFC为强迫对流无量纲数。
由于重力自循环系统涉及到气液两相流动，物理机制较复杂，因而国内外的研究大多用CFD软件单独模拟热沉部分的流动，计算耗时较长，且不能得到系统的整体特性。
本文提出了针对重力自循环系统的一维仿真模型，研究系统在不同加热功率下的总体流量及压降特性。作者研究了不同系统结构参数和边界条件下系统的性能，可为一般重力自循环系统的设计提供参考。
1 数学模型及边界条件
由于计算工况中干度较小[见图3(c)]，因而采取一维均相流动模型，将气液混合物看作拥有平均物理性质的赝流体。这里进行简化，认为沸腾起始点(ONB)点之前，流体干度x等于0，ONB点后流体处于饱和状态。
ONB点位置与壁温密切相关，因而需要选取合适的换热系数得到壁温，采用式(2)计算[11]
(2)
式中，Pr为流体的普朗特数。则壁温可计算如下
(3)
其中，l为管径d。得到壁温后，即可根据式(4)确定ONB点位置
(4)
流动满足质量守恒方程
(5)
(6)
其中：ρm、ρl、ρg分别为混合密度、液氮密度、气氮密度；A为管截面积；v为流体速度。
能量守恒方程
(7)
其中：q为热流密度；D为管径；G为系统流量；H为流体的比焓；v为流体的速度；g为重力加速度；z为竖直高度；θ为管倾斜角
H=cpl·Tf·(1-x)+r·x
(8)
式中：cpl为流体比热容；Tf为流体温度；r为汽化潜热；x为干度。起泡前，基本可以认为x等于0，起泡后，由于流体已达饱和状态，可以认为液体部分温度Tf=Tsat。根据能量守恒，干度计算如下
(9)
其中，zFDB表示起泡点所在高度。
流体温度可由式(7)推出
(10)
流动达到稳定时，系统的重力压差与流动阻力达到动态平衡，因而正确计算系统各部分的压降至关重要。系统的阻力包括两部分：沿途阻力损失与局部阻力损失。
沿途压降可根据动量守恒方程计算
(11)
两相流的压降可分为3部分：重力压降、摩擦压降与加速压降。压降根据Martinelli-Chisholm关联式[12]计算，邓东等证实此关联式适用于液氮[13]。即：
(12)
(13)
(14)
(15)
其中：μl、μg分别为液氮与气氮的动力黏度；ρl、ρg分别为液氮与气氮的密度；x为干度。
重力压降与加速压降计算如下：
(16)
(17)
其中：G为系统流量；A为管道截面积。
局部压力损失系数为：
(18)
系统在流动达到稳定时，整体的动力压头与总阻力是相等的，因而出口处流体的静压强应等于外界环境压强，否则系统流量仍会波动，如图1所示，即有：
pout=pa
(19)
计算中，出口处压强pout根据入口压强与沿途压降相加得到，环境压强pa是已知的。
式(3)、式(5)、式(9)、式(10)、式(11)、式(19)六式组成的方程组中，在z
本研究采用Baudouy的实验结果[10]验证模型的正确性，相关边界条件如下：①pa取1个大气压；② 热流密度q变化范围为2～12 kW/m2；③ 管长、管径等结构参数见第2节。
2 实验装置介绍与算法说明
本研究首先采用Baudouy的实验数据验证模型的正确性，其实验装置工作状态如图1所示，实验装置由循环回路与外部的绝热罩两部分组成。循环回路的上升段与下降段通过U型与液氮储箱连接，下降段管径为40 mm，上升段管径为10 mm。液氮储箱直径为0.45 m，高0.3 m。上升段长1.45 m，其中受热部分位于上升段的下部，长0.95 m，非受热段长0.5 m。实验通过安置文氏流量计测量系统流量，文氏流量计长0.4 m，入口直径为40 mm，喉部直径为10 mm。实验装置的具体结构参数见文献[10,14]。
计算时，初设流量后进行迭代对数学模型进行求解，其流程如图2所示，算法在Matlab中实现。
3 结果分析
3.1 模型正确性验证
以Baudouy的实验系统为对象进行仿真计算，计算值与实验值的对比如图3所示。
由图3可见，随着热流密度增大，计算流量先增大，达到极大值后开始减小，流量的相对误差绝大部分均在±10%以内；受热段总压降随着热流密度的升高不断减小，压降的相对误差绝大部分在±10%以内。流量与压降随热流密度的变化趋势均与实验结果相符，且相对误差较小，验证了仿真模型的正确性。
3.2 管径对系统工作特性的影响
热沉管径D对流体的流动与换热特性有着重要影响。在其他参数不变时，研究不同管径对系统工作特性的影响。
图4(a)中，不同管径下，流量随热流密度升高的变化趋势是一致的，都是先增大后减小。在同一热流密度下，大管径对应着较大的流量。而其他参数相同时，大管径所对应的单位长度上的换热量较大，因而对应着较大的干度，所以重力压差较大；另一方面，摩擦压降与管径成反比关系，两种作用使得大管径系统流量更大。
随着管径增大，使系统达到最大流量的临界热流密度也在逐渐增加大。在管径为0.005 m时，临界热流密度在 2 000 W/m2左右，而在管径为0.02 m时，临界热流密度已经超过了10 000 W/m2。图4(c)展示了不同管径下系统出口处的干度值，在管径较小时，受热管内干度较大，摩擦阻力较大，使得流量在较小的热流密度下达到极大值。
3.3 受热段长度对系统工作特性的影响
液氮在受热段吸热汽化，受热段长度决定着液氮总吸热量，进而影响系统内气相组分大小。在Baudouy的实验装置基础上，保持上升段长度为1.45 m不变，取受热段起始点高度H0=0，改变受热部分的长度，计算结果如图5所示。
对比图5发现，受热长度对系统流量的影响与分受热长度大小有关。在受热长度较小(Lh/L<0.4)时，受热长度增加使系统流量增大；在受热长度满足0.6>Lh/L>0.4时，受热长度增加在低热流密度下使系统流量增大，在高热流密度下使系统流量减小。在受热长度较大(Lh/L>0.4)时，受热长度增加反而使系统流量减小。
上升段内的平均干度受热流密度与加热长度的共同影响。存在一个临界的平均干度xc，在平均干度小于xc时，上升段内重力压降起主要作用，在平均干度大于xc时，摩擦压降替代重力压降起主要作用。在受热长度较小时，上升段内干度较小，因而受热长度增加使平均干度增加，上升段与下降段的重力压差随之增大，导致系统流量增大。在受热长度较大时，上升段内平均干度较大，而摩擦压降随干度的增加而增加，因而系统流量反而随受热长度的增大而减小。在受热长度处于两者之间时，受热长度对流量的影响在高热流密度区与低热流密度区截然相反。
3.4 受热段位置对系统工作特性的影响
加热段长度一定时，加热位置不同会改变流动的阻力特性。在Baudouy的实验装置基础上，保持上升段长度为1.45 m不变，取Lh/L=0.2，改变受热段起始点高度H0，结果如图6所示。
总体上，受热位置上升会使得系统流量减小。分析认为，系统循环流动的动力在于流体受热汽化后在下降段与上升段之间产生了密度差，密度差导致重力压差。而两相区的长度越大，重力压差越大。因而，加热段位于下部时，液氮汽化后流经的长度更大，使得上升段与下降段的压强差更大，流量更大。图6(b)证实了这一点。
干度与受热长度近似成正比关系，只有受热长度达到一定值时，干度才能产生较大影响。因此可以推断：受热长度越大，受热位置上升使系统流量减小的趋势越明显。
3.5 上升段总管长对系统工作特性的影响
在Baudouy实验装置中，非受热段处于受热段的上方，因而上升段与下降段的重力压差由受热段与非受热段共同提供，因此上升段的总管长对系统工作特性有着重要影响。Baudouy的实验装置中，系统上升段总长度为1.45 m，在实际工程中管长有时候会达到数十米乃至百米[7]。为研究管长对系统工作特性的影响，这里保持系统的受热段长度为0.95 m不变，将上升段的管长Lup增大，同时相应地增大下降段的管长Ldown，计算结果如图7所示。
总管长的增加对系统流量的影响需要分低热流密度区与高热流密度区讨论。热流密度较低时，管长增加对系统流量的影响不大，热流密度较高时，管长增加会使得系统流量有一定程度的增大。分析认为，在热流密度较低时，上升段内主要是单相流，这使得上升段与下降段间压差减小，使得系统在较小的流量下即可实现换热和动量平衡。在热流密度达到3 000 W/m2以上时，由于受热管内干度与热流密度近似成正比关系，受热段内干度足够大，使得上升段内重力压降起主要作用，因而管长增加会使得系统流量增大。
4 结论
1) 仿真模型的仿真结果与Baudouy的实验数据相当一致，说明该模型能对一般重力自循环系统的稳态特性进行正确计算。
2) 仿真研究发现，系统的结构参数如管径等对系统的流量与压降特性有重要影响。
3) 本文研究工作可用于估算已有的重力自循环热沉系统的性能，也可用于设计新系统。
参考文献：
[1] 单巍巍，刘波涛，丁文静,等.重力式自循环系统中热沉结构设计方法研究[J].航天器环境工程,2010,27：489-492.
[2] DMITRIY PODKORYTOV,LEONID TIMKIN,VLADIMIR CHECHOVICH.The investigation of flow characteristics and heat transfer at the nitrogen natural circulation in capillary[J].Cryogenics,1992,32(1):316-319.
[3] YEON SUK CHOI,HO-MYUNG CHANG,STEVEN W.Van Sciver.Performance of extended surface from a cryocooler for subcooling liquid nitrogen by natural convection[J].Cryogenics,2006,46(5):396-402.
[4] KIM M J,CHANG H.M.Natural circulation loop of subcooled liquid nitrogen[C]//Advances in Cryogenic Engineering.Chattanooga,2008.
[5] CHANG HO-MYUNG,KIM MIN-JEE,SIM JUNG WOOK.A compact cryocooling system with subcooled liquid nitrogen for small HTS magnets[C]//20th International Conference on Magnet Technology.Philadelphia,2007.
[6] YU.IVANOV,A.RADOVINSKY,A.ZHUKOVSKY,et al.Compact counter-flow cooling system with subcooled gravity-fed circulating liquid nitrogen[J].Physica C:Superconductivity,2010,470 (20):1895-18989.
[7] IVANOV Y,RADOVINSKY A,ZHUKOVSKY A,et al.A compact cooling system for HTS power cable based on thermal siphon for circulation of LN2[C]//Transactions of the Cryogenic Engineering Conference—cec:Advances in Cryogenic Engineering.American Institute of Physics,2010:865-870.
[8] YAMAGUCHI S,IVANOV Y,SUN J,et al.Experiment of the 200-Meter Superconducting DC Transmission Power Cable in Chubu University[J].Physics Procedia,2012,36:1131-1136.
[9] YURY IVANOV,HIROFUMI WATANABE,MAKOTO HAMABE,et al.Observation of the thermosiphon effect in the circulation of liquid nitrogen in HTS cable cooling system[J].Physics Procedia,2012,27:368-371.
[10] BAUDOUY.Experimental study of a nitrogen natural circulation loop at low heat flux[C]//Joint Cryogenic Engineering Conference.Tucson,2009.
[11] 李友荣.高等传热学[M].北京：科学出版社，2013.
[12] 徐济鋆.沸腾传热和气液两相流[M].北京:原子能出版社，2001.
[13] DENG D,XIE S W,WANG R S.Two-phase flow pressure drop of liquid nitrogen boiling in the straight section downstream of U-bend[J].Journal of Shanghai Jiaotong University (Science),2014,19(4):495-50187.
[14] BENKHEIRA L,BAUDOUY B,SOUHAR M.Heat transfer characteristics of two-phase He I (4.2 K) thermosiphon flow[J].International Journal of Heat and Mass Transfer,2007,50(17):3534-3544.
