热负荷变化对加热管内正十烷的振荡数值研究

田 静 ，侯 宇 ，周伟星，肖雪峰 ，王广飞

（1．中国民航大学航空工程学院，天津 300300；2.哈尔滨工业大学a.基础与交叉科学研究院；b.能源科学与工程学院，哈尔滨 150006）

碳氢燃料的再生冷却技术是超燃冲压发动机最为重要和有效的保护措施[1-2]，碳氢燃料在进入燃烧室之前要通过流动通道对发动机壁面进行吸热冷却，这对提高发动机的性能和燃料的充分利用十分有利[3]。在碳氢燃料吸收热量升温的过程中温度升高超过临界温度，并且温度继续升高出现化学裂解反应，进一步吸收热量使得发动机壁面温度进一步下降[4-5]。由于燃料进入裂解区，物理性质和化学性质发生十分剧烈的变化，以往的研究大都集中于燃料在管道内流动换热的稳态研究，本文只针对正十烷非裂解状态条件下展开动态特性研究。

李勋峰等[6]采用数值模拟的方法对超临界压力下燃油管道的传热特性进行研究，发现煤油的传热与近壁面区域的流动状态有关。赵国柱等[7-8]通过数值计算与实验研究了正十烷在超临界压力下的流动传热，发现目前的传热公式在正十烷的临界区域与计算结果相差很大，并且进一步通过实验研究了正十烷在超临界压力条件下在湍流区的传热采用多元线性回归方法，得出了正十烷在湍流区的传热关联式并且验证其适用性。姜蕾等[9]对碳氢燃料进行了变压力过程中传热现象的研究，发现变压力过程中壁面温度的非单调特性，并且在超临界压力下动态变化时传热弱于稳态传热。Wang等[10]和Hua等[11]分别对低温甲烷和正庚烷在水平管内的超临界湍流传热进行了分析，研究表明：高压力时二者的换热性能良好；在接近临界压力的情况下，当壁温处于拟临界温度附近时，出现了传热恶化现象。Zhou等[12]对超临界戊烷的对流换热研究中也发现了相似的规律，并给出了传热恶化发生的临界热流密度关系式。Brad Hitch[13]研究发现超临界条件下，Methyhyclohexane（MCH）在不同压力和温度区间出现了高频和低频两种振荡现象，并提出在管内增加扰流器是一种有效的抑制方法。于彬等[14]建立了0维平衡流模型，对不同裂解特性燃料进行研究，研究表明：供油系统存在两个不稳定区即临界区和裂解区，燃料密度和裂解速率是不稳定现象的主要原因。

本文采用了FLUENT计算平台对N-S（Navier-Stokes）方程进行离散求解，且正十烷的物性计算与入口质量流量和壁面热负荷变化进行耦合计算，湍流模型采用了k-ω SST双方程模型[15]。为验证计算方法的准确性，进行了网格无关性验证并且与文献实验数据进行了对照，采用验证后的计算方法对壁面热负荷动态变化进行了数值计算。对处于拟临界温度范围内的正十烷在壁面热负荷变化条件下的不稳定现象进行了研究。

1 计算方法

1.1 数学模型

加热管道采用高温合金钢（如图1所示），其密度ρ=8 810 g/m3，定压比热容 Cp=142.202 6 J/kg·K，导热系数k=16.75 W/m·K，壁厚d=0.5 mm，直径D=2 mm，入口绝热段l1=100 mm，使得管内湍流充分发展，加热段l2=1 000 mm，出口绝热段l3=100 mm，加热段管道壁面采用加载轴对称恒定体热流密度。假设加热段没有向外界散失热量，管道壁面传入的热量都被流体完全吸收，忽略管道的轴向导热。根据入口边界条件，管道内的G（r格拉晓夫参数）同Re数的平方之比的量级比大约在10-5～10-4，表明重力对换热影响很小。

本文算例参数设置如下：

1）入口采用质量流量边界条件，正十烷的入口初始质量流量为0.004 9 kg/s；

图1 加热管道物理模型Fig.1 Physical model of heating pipe

3）出口采用压力出口边界条件，压力为3 MPa；

4）设置壁面热负荷边界条件。

通过大量实验可以发现：壁面热负荷变化基本有两种，一种是近似阶跃变化，另一种则是近似线性变化。因此本文采用这两种热负荷变化过程作为仿真对象。如图2所示，其中a代表线性热负荷变化（line），热负荷随时间变化关系为W=1 624+150.6 t，0≤t≤1.5 s；b代表热负荷阶跃型（step），热负荷阶跃变化如表1所示。

图2 热负荷变化过程Fig.2 Heat load change process

表1 热负荷阶跃型参数Tab.1 Heat load step parameters

1.2 数值方法

本文采用有限体积法对N-S方程进行离散求解，采用SIMPLE算法作为压力-速度耦合方法求解，入口雷诺数大于圆管湍流最低值2 300并且经过入口段充分发展，流体基本进入湍流旺盛区，故湍流模型采用k-ω SST模型进行传热模拟有较好的效果。

在直角坐标系中，以张量形式表示的可压缩流体的时均控制微分方程如下：

1）连续方程

2）动量方程

3）能量方程

4）湍动能输运方程

标准模型的k-ω SST湍动输运方程为

5）比耗散输运方程

在上述方程中，Pk、Dk分别为k的生成项和耗散项；Fl为混合函数；αk、αω分别为k和 ω 的等效耗散系数。由于是圆管模型则以管道的中轴线为对称轴，计算采用二维轴对称简化模型，减少计算量。

1.3 物性计算

由于正十烷在管道内的流动状态处于超临界压力下，其物性受温度变化的影响十分明显，尤其是在临界点的附近，其可压缩性较强，因此采用PR状态方程[16]确定体积、压力、温度的变化，定压比热容采用多项式拟合的方式。为了准确模拟正十烷物性变化，粘度和导热系数均采用Chung等[17-18]的计算方法。物性计算方法通过编写程序，后通过导入FLUENT计算软件user define function，用以研究正十烷在加热管道的动态特性。

1.4 网格无关性验证

由于超临界条件下正十烷的物性受温度变化的影响较大，为了比较准确的计算传热采用k-ω SST模型，要求网格保证在粘性底层有足够的点，因此计算网格y+＜5[5]的大小与网格质量有关。如表2所示，对4组不同网格进行验证。

表2 网格无关性验证Tab.2 Grid independence verification

可发现4种不同网格计算所得的壁面温度相差非常小，约为0.2%，为节省计算资源，故选标号4作为本文所使用的网格，网格数为60×2 200。

1.5 计算验证

为计算非稳态条件下的数值模拟，首先验证了稳态时模型的准确性和数值方法的准确性。为此与文献[7-8]的稳态实验数据进行对照验证，如图3所示。

图3 文献验证Fig.3 Document verification

文献[19]按照雷诺数划分湍流模型，当Re＜2 300时为层流区，2 300≤Re≤10 000时为过渡区，Re＞10 000时为湍流旺盛区。由于计算采用k-ω SST湍流模型，根据文献[7]其管长l≤0.4 m，在实验中管道内正十烷尚处于层流区和过渡区，故传热较慢，壁面温度较高。而当l≥0.4 m时，计算壁面温度与实验数据有较好的吻合。此计算模型与数值计算方法适合高雷诺数条件下正十烷的流动换热模拟。

由于超临界状态下正十烷的物性变化受温度影响较大，为了说明本文所采用计算模型的准确性，本文计算了3 MPa压力条件下热负荷a调节方式，变化初始时刻轴线处正十烷的Re变化曲线如图4所示，正十烷线性热负荷变化初始时刻，当l≥0.2 m时，管内正十烷都进入湍流旺盛区，随着壁面热负荷的继续增大，管内正十烷吸热升温Re持续增大，因而管内绝大部分正十烷都进入湍流旺盛区，故数值计算采用的湍流模型能够较好地用于研究正十烷的动态特性变化过程。

图4 加热管内Re沿着管长变化Fig.4 Re change along tube length in heating tube

2 数值模拟结果及其分析

2.1 热负荷变化速率对出口参数变化的影响

壁面热负荷按照线性a与阶跃b1方式变化（如图2（a）所示），出口质量流量变化曲线如图5所示。

图5 出口质量流量变化比较Fig.5 Comparison of exit quality and flow

热负荷按照阶跃型与线性变化都引起进出口质量流量不守恒现象。如图5所示，阶跃型热负荷调节方式（b1）出口质量流量超调现象比较显著，在较短时间内出口质量流量大幅高频振荡，而后迅速下降趋平稳，而线性热负荷调节方式（a）出口质量流量呈现一个“梯形”变化方式，出口流量超调量明显下降；壁面热负荷变化越快，其出口质量流量超调值越大。

2.2 热负荷阶跃小幅值对流体温度场的影响

为了研究流体跨越拟临界温度区的振荡现象，热负荷变化过程为b2，如图2（b）所示。首先研究了非拟临界区正十烷的温度变化规律，如图1所示取管道中轴线两点 c（0.815 m，0）、d（0.965 m，0），如图 6 所示：c点温度变化规律，当t≤0.005 s时其温度持续降低；0.005 s＜t≤0.008 s时温度上升非常缓慢，而后当t＞0.008 s时温度上升速率增加，没有出现反复振荡现象。

图6 热负荷阶跃对非拟临界区正十烷温度变化的影响曲线Fig.6 Effect curve of heat load step on n-decane temperature in non-quasi-critical region

图7为同一热负荷变化过程b2，d点处于拟临界温度范围，研究其温度的变化规律，其在t=0.001～0.01 s内温度下降2.4 K，而后在t=0.011～0.021 s温度上升约1.9 K，而当t=0.022～0.025 s下降约0.3 K，其后温度持续上升。对比图6和图7可以发现，对于同一热负荷小幅值阶跃变化过程中，处于拟临界范围的正十烷，其温度的振幅和振荡的频率都比非拟临界区域的正十烷高。

图7 热负荷阶跃对拟临界区正十烷温度变化的影响曲线Fig.7 Effect curve of heat radiation step on n-decane temperature in quasi-critical region

2.3 热负荷阶跃大幅值对流体温度场的影响

如图1所示取管道中轴线上两点c（0.815 m，0）、d（0.965 m，0），如图8所示：c点的温度处于非拟临界区，对热负荷阶跃变化过程b2和b3进行比较，c点温度变化规律在b2、b3变化过程中基本上是一致的，温度先下降而后上升。当t≤0.017 s时，b2过程c点温度高于b3过程。而当t＞0.017 s后，b3过程c点温度逐渐超过b2过程，c点温度增加速率b3过程高于b2过程。

图8 热负荷幅值阶跃大小对非拟临界区正十烷温度变化影响Fig.8 Effect of thermal load amplitude step size on n-decane temperature in non-quasi-critical region

图9表示热负荷阶跃变化过程b2、b3，d点处于拟临界温度范围，研究其温度变化规律。由于物性变化十分剧烈，可以发现在b2、b3变化过程中该点温度有2次振荡。t=0～0.01 s时，温度均先下降且下降幅值随着阶跃幅值的增大而增大；当t=0.01～0.02 s时温度都回升，此过程末热负荷阶跃幅值大的温度略低于阶跃幅值小的变化。

图9 热负荷幅值阶跃大小对拟临界区正十烷温度变化影响Fig.9 Effect of thermal load amplitude step size on n-decane temperature in quasi-critical region

阶跃为b2时，其温度在t=0.022～0.025 s时再次降低，但其温度下降幅值较初次低；当t=0.025 s后，温度又单调增加。阶跃为b3时，在t=0.022～0.03 s时，温度降幅较大；而t=0.03 s后，温度单调增加；当t≤0.045 s时，热负荷幅值小阶跃条件下d点的温度始终高于热辐荷阶跃幅值较大的变化情况。

2.4 热负荷阶跃幅值变化对流体速度场影响

热负荷幅值变化大小对出口正十烷的速度影响如图10所示，其中选取出口截面两点分别为e（1.2 m，0.008 m）和（f1.2 m，0.002 m）。对比b2和b3过程可看出，同一热负荷变化过程f点的平均速度要高于e点正十烷的平均速度；流体的速度沿着半径指向轴心不断增加；热负荷阶跃幅值越大，f点处的平均速度的振幅越大；f和e点速度的平均值随时间变化的规律基本一致，均是先增加，而后下降，再增加到低幅振动区，随后继续增加，即热负荷阶跃幅值越大时，出口流体的速度振荡越大。

图10 热负荷幅值变化大小对出口正十烷的速度变化影响Fig.10 Effect of heat load amplitude change on n-decane at exit

3 结语

由于正十烷的温度升高物性发生剧烈变化，其可压缩性显著提高，从而导致了燃油不稳定现象。

1）非拟临界温度区域内的正十烷在热负荷阶跃变化时其温度变化为先下降而后上升，且温度先下降的趋势沿着加热管道逐步增强；阶跃幅值越大对处于非拟临界温度区的正十烷的温度影响越大。

2）热负荷变化快慢对处于拟临界温度区的正十烷的振荡是显著的，热负荷增加速率越快，振荡现象越显著；壁面热负荷变化越快，其出口质量流量超调值越大。

3）热负荷阶跃幅值越大，拟临界温度范围内正十烷的温度振荡现象越显著；正十烷处于拟临界温度范围内，热负荷阶跃幅值变化较大的过程，其温度在振荡期内始终低于热负荷阶跃幅值较小的变化过程。

4）热负荷阶跃值越大，速度场的振幅越大，频率也越高。

本文对超临界压力条件下燃油供给过程中不稳定现象的成因进行了研究，提供了一定参考价值。

[1]HUANG H,SPADACCINI L J,SOBEL D R.Fuel-cooled thermal management for advanced aero-engines[J].J Eng Gas Turbines Power,2004,126:284-293.

[2]SOBEL D R,SPADACCINI L J.Hydrocarbon fuel cooling technologies for advanced propulsion[J].J Eng Gas Turbines Power,1997,119:344-351.

[3]LINNEDL,MEYERML,EDWARDST,etal.Evaluationofsupercritical fuel injection on cycle performance of pulsed detonation engine[J].J Propulsion Power,2007,23(4):748-755.

[4]HELFRICH T M,KING P I,HOKE J L,et al.Effect of supercritical fuel injection on cycle performance of pulsed detonation engine[J].Propul－sion Power,2007,23(4):748-755.

[5]LI X,HUAI X,CAI J,et al.Convective heat transfer characteristics of China RP-3 aviation kerosene at supercritical pressure[J].Appl Therm Eng,2011,31(14):2360-2366.

[6]李勋峰,仲峰泉,范学军,等.超临界压力条件下航空煤油圆管流动和传热的数值研究[J].推进技术,2010,31:467-472.

[7]赵国柱,宋文艳,张若凌,等.超临界压力下正十烷流动传热的数值模拟[J].推进技术,2014,35(4):537-543.

[8]张 磊,乐嘉陵,张若凌,等.超临界压力下湍流区碳氢燃料传热研究[J].推进技术,2013,34(2):225-229.

[9]姜 蕾,刘朝辉,毕勤成,等.碳氢燃料变压力过程中传热研究[J].推进技术,2014,35(7):966-971.

[10]WANG Y Z,HUA Y X,MENG H.Numerical studies of supercritical turbulent convective heal transfer of cryogenic propellant methane[J].Journal of Thermo Physics and Heat Transfer,2010,24(3):490-500.

[11]HUA Y X,WANG Y Z,MENG H.A numerical study of supercritical forced convective heal transfer of n-Hep-tane inside a horizontal miniaturetube[J].JournalofThermalSupercritical Fluids,2010,52(1):36-46.

[12]ZHOU W X,BAO W,QIN J.Deterioration in heal transfer of endothermal hydrocarbon fuel[J].Journal of Thermal Science,2011,20(2):173-180.

[13]BRAD H.Enhancement of Heat Transfer and Elimination of Flow Oscillations in Supercritical Fuels,AIAA-98-3759[R].

[14]于 彬,周伟星,于文力,等.燃料裂解特性对供油系统稳定性的影响[J].推进技术,2013,34(12):1703-1707.

[15]MENTER F R.Two-equation eddy-viscosity turbulence models for engineering applications[J].AIAA Journal,1994,32(8):1598-1605.

[16]冯 宇.超临界碳氢燃料流动裂解藕合及传热特性的数值研究[D].哈尔滨:哈尔滨工业大学,2014.

[17]CHUNG T H,LEE L,STARLING K E.Applications of kinetic gas theories and multiparameter correlation for prediction of dilute Gas viscosity anthermal conductivity[J].Industrial&Engineering Chemistry Fundamentals,1984,23(1):8-13.

[18]CHUNG T H,AJLAN M,LEE L L,et al.Generalized multiparameter correlation for nonpolar and polar fluid transport properties[J].Industrial&Engineering Chemistry Fundamentals,1988,27(4):671-679.

[19]陶文栓.数值传热学[M].2版.西安:西安交通大学,2006.