查询字段 检索词
  推进技术  2018, Vol. 39 Issue (9): 2028-2034  DOI: 10.13675/j.cnki.tjjs.2018.09.013
0

引用本文  

初敏, 陈思员, 史可天, 等. 可压缩算法在燃油冷却数值研究中的应用[J]. 推进技术, 2018, 39(9): 2028-2034.
CHU Min, CHEN Si-yuan, SHI Ke-tian, et al. Application of Compressible Algorithm in Numerical Simulation of Fuel-Active-Cooling[J]. Journal of Propulsion Technology, 2018, 39(9): 2028-2034.

通讯作者

艾邦成,男,博士,研究员,研究领域为高超声速飞行器热防护设计。E-mail:aimen011@126.com

作者简介

初敏,男,博士,研究领域为高超声速飞行器热防护设计。E-mail:chumin_caaa@126.com

文章历史

收稿日期:2017-09-08
修订日期:2017-11-02
可压缩算法在燃油冷却数值研究中的应用
初敏 , 陈思员 , 史可天 , 艾邦成     
中国航天空气动力技术研究院,北京 100074
摘要:为深入研究超燃冲压发动机的燃油主动冷却过程,在可压缩算法框架下,引入预处理技术、低耗散数值离散格式、一般流体状态方程、高压物性模型、燃油裂解模型等,发展了一套适用于求解一般流体低速可压缩问题的数值模拟方法。利用电加热管试验对数值方法进行了验证,并定量研究了PR-3燃油的冷却需求量及平板面板冷却槽道布置的影响。结果表明,数值结果与试验符合很好,证明当前数值方法在燃油冷却数值模拟中是准确有效的。针对350mm×90mm×6mm的平板面板冷却过程,需要约10g/s的燃油可确保1WM/m2热流下的防热安全。槽道布置对温度场影响明显,两侧进油方式比单侧进油的最大温差低约150K,均温性更好,更有利于降低结构热应力。
关键词燃油主动冷却    数值模拟    可压缩算法    预处理    
Application of Compressible Algorithm in Numerical Simulation of Fuel-Active-Cooling
CHU Min, CHEN Si-yuan, SHI Ke-tian, AI Bang-cheng     
China Academy of Aerospace Aerodynamics, Beijing 100074, China
Abstract: In order to study the fuel-active-cooling of scramjet, a numerical simulation method suitable to solve the low-speed compressible problems of general fluid was developed. Based on the framework of compressible algorithm, the method introduced precondition technology, low dissipation numerical scheme, general fluid state equation, physical property model of high pressure and fuel pyrolysis model. The results of electrically heated tube test were adopted to verify the physical models and algorithms. The mass flow rate of RP-3 fuel for cooling process of flat panel was analyzed quantitatively and the effects of the cooling channel arrangement on cooling process were studied. The results show that the physical models and algorithms are effective on fuel cooling. The mass flow rate of RP-3 should be more than 10g/s to ensure the thermal protection safety under 1MW/m2 heat flux when the size of flat panel is 350mm×90mm×6mm. The channel arrangement has obvious influence on the temperature field. The maximum temperature difference when there are two oil inlets is 150K lower than that when there is only one oil inlet. Better isothermal performance means it is beneficial to the decrease of structure thermal stress.
Key words: Fuel-active-cooling    Numerical simulation    Compressible algorithm    Precondition    
1 引言

作为吸气式巡航高超声速飞行器的主要动力装置,超燃冲压发动机的热防护是制约飞行器发展的一项关键技术。为了满足长时间飞行和多次可重复使用的需求,具有非烧蚀特点的主动冷却方式已逐渐成为发动机热防护的优选方案。

采用燃油作为冷却剂对发动机进行主动冷却的方式,在液体火箭发动机领域较为成熟,由于携带流量大,燃油在冷却过程中的温升相对较小,其流动传热过程相对简单。而高超声速飞行器为了减少质量惩罚,只携带有限量的燃油,即使采用热沉更大的吸热型碳氢燃料,如美国产JP-7和中国产RP-3[1]等,仍需充分利用燃油的热沉,导致燃油在冷却过程中的温升较大,超过燃油的裂解温度,从而发生热裂解及结焦等一系列现象。此外,为了避免温度过高发生沸腾,冷却通道内的燃油一般处于超临界压力下。因此,超燃冲压发动机的燃油主动冷却过程,实际上是燃油由过压液态逐渐吸热升温变为超临界态,并进一步裂解为小分子物质的过程,其中还可能伴随着燃油的氧化结焦和裂解结焦,流动传热过程非常复杂[2]。目前试验手段难以对通道内部进行细致观察,需借助于计算流体力学方法开展更加深入的研究。

国内在针对燃油主动冷却问题开展数值研究时,较多采用了SIMPLE系列的算法。清华大学的王夕[2]采用SIMPLEC算法研究了冷却通道入口汇流槽内的流动结构对RP-3煤油传热过程的影响,指出汇流槽内存在非定常流动过程,会导致各通道内的流量差异达到10%。西北工业大学的赵国柱等[3]采用SIMPLE算法研究了正癸烷的超临界传热过程,计算发现常见的对流传热计算公式在临近点附近误差明显。赵国柱等还针对RP-3煤油[4]的流动传热过程开展了研究,着重分析了二次反应[5]对裂解的影响。中国空气动力研究与发展中心的张若凌等[6]采用SIMPLE算法结合煤油多步裂解反应机理开展了燃油主动冷却过程的设计分析,计算结果与试验结果符合较好。天津大学的张强强[7]采用SIMPLEC算法研究了正癸烷在超临界条件下的传热特性,通过引入大分子链烯烃的二次反应,提出了改进的正癸烷裂解模型,用于模拟正癸烷在高裂解率下的反应过程,效果较好。文献[8, 9]采用SIMPLEC算法研究了冷却通道横截面参数非均匀性及通道纵横比对燃油冷却过程的影响,指出增大纵横比可以增强传热,但过大的纵横比会导致热分层,使煤油热沉利用不充分。

超燃冲压发动机的燃油主动冷却过程是一个低速可压缩的特殊过程,燃油升温后会进入超临界态并发生裂解,组分、动量、能量方程强烈耦合,表现出较强的压缩性。尽管以SIMPLE为代表的不可压算法已推广至可压缩问题的求解,但在压缩性较强的问题中,可压缩算法具有更好的收敛性。另外,目前的不可压缩算法和可压缩算法是相互独立的,而用于高超声速数值模拟的主要是基于密度求解的可压缩算法。在可压缩算法框架下,将其扩展至低速不可压问题的求解,有助于实现算法的统一,便于开展多物理过程的耦合计算。因此,本文对可压缩算法进行扩展,将其用于燃油主动冷却过程的数值模拟,对可压缩算法在此问题下的有效性进行了验证,并进一步对平板面板的冷却过程进行了分析。

2 物理模型和计算方法 2.1 控制方程及数值离散方法

流体控制方程采用的是多组分形式下的NS方程,其张量形式如下

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {\rho {Y_j}} \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho {Y_j}{u_i}} \right) = \frac{\partial }{{\partial {x_i}}}\left( {{D_{j{\rm{m}}}}\rho \frac{{\partial {Y_j}}}{{\partial {x_i}}}} \right) + {{\dot \omega }_j}, \\ j = 1, {\rm{}}2, \cdots , ns\\ \frac{\partial }{{\partial t}}\left( {\rho {u_j}} \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho {u_j}{u_i} + {\delta _{ij}}p} \right) = \frac{{\partial {\tau _{ij}}}}{{\partial {x_i}}}\\ \frac{\partial }{{\partial t}}\left( {\rho E} \right) + \frac{\partial }{{\partial {x_i}}}\left[ {\left( {\rho E + p} \right){u_i}} \right] = \\ \frac{\partial }{{\partial {x_i}}}\left[ {{\tau _{ij}}{u_j} - {q_i} + \rho \sum\limits_{j = 1}^{ns} {D_{j{\rm{m}}}}{h_j}\frac{{\partial {Y_j}}}{{\partial {x_i}}}} \right] \end{array} $ (1)

式中ns为组分个数,ρ为流体密度,Yj为第j种组分的质量分数,uj为流体的速度分量,p为压力,E为流体总能量,H为总焓,τij为应力张量,Djm为第j种组分的平均分子扩散系数,qi为热传导项,$ {\dot \omega _j}$为化学反应源项。流动状态对传热过程影响很大,采用剪切应力输运(Shear Stress Transport,SST k-ω)模型[10]模拟湍流。

控制方程中的时间项采用LUSGS方法[11]进行离散。对于空间项,由于求解的是一般流体在低速下的流动传热问题,可压缩算法中经典的迎风系列格式并不能直接使用,采用Edwards提出的针对一般流体在全流速下的低耗散格式LDFSS(Low-Diffusion Flux Splitting Scheme)[12]进行离散。为进一步加速收敛,引入了多重网格技术。

2.2 一般流体状态方程

燃油在冷却通道内会经历亚临界到超临界态的转变,流体性质从液态向气态过渡,理想气体状态方程不再适用,需要引入可以同时描述气态和液态流体热力学关系的状态方程。立方型状态方程由于物理意义明确、形式简单,便于推导各种热力学关系式,应用较多。在常见的几种立方型状态方程中,Soave-Redlich-Kwong(SRK)方程对小分子碳氢化合物较为准确,而Peng-Robinson(PR)方程对大分子碳氢化合物较为准确[13, 14]。在燃油冷却过程中,燃油会发生裂解,由初始的大分子物质生成氢气、乙烯等一系列小分子物质,单一的状态方程难以对整个过程进行准确的描述。Cismondi等[14]通过引入系数δ1δ2,提出了RK-PR状态方程,将SRK方程和PR方程结合起来,从而提高了状态方程的适用范围。RK-PR状态方程的形式如下

$ p = \frac{{RT}}{{V - b}} - \frac{{a\left( T \right)}}{{\left( {V + {\delta _1}b} \right)\left( {V + {\delta _2}b} \right)}} $ (2)
$ \left\{ {\begin{array}{*{20}{l}} { - {\delta _1}{\delta _2} = c}\\ {{\delta _1} + {\delta _2} - 1 = c} \end{array}} \right. $ (3)

c取0时,RK-PR方程退化为SRK方程,取1时,则退化为PR方程。其它参数形式可参考原文献。当流体为混合物时,采用如下混合规则

$ \begin{array}{*{20}{l}} {a = \sum\limits_{i = 1}^n \sum\limits_{j = 1}^n {X_i}{X_j}{a_{ij}}\;\;\;\;b = \sum\limits_{i = 1}^n {X_i}{b_i}}\\ {{\delta _1} = \sum\limits_{i = 1}^n {X_i}{\delta _{1, i}}\;\;\;\;\;\;\;\;\;\;{\delta _2} = \sum\limits_{i = 1}^n {X_i}{\delta _{2, i}}} \end{array} $ (4)

式中X为摩尔浓度。根据状态方程形式,可推导压力与温度、密度、组分质量分数的热力学关系式。

$ {\left( {\frac{{\partial p}}{{\partial T}}} \right)_{{\rho _i}}} = \frac{{R\rho }}{{M - b\rho }} - \frac{{{\rho ^2}}}{{\left( {M + {\delta _1}b\rho } \right)\left( {M + {\delta _2}b\rho } \right)}}\frac{{\partial a\left( T \right)}}{{\partial T}} $ (5)
$ {\left( {\frac{{\partial p}}{{\partial \rho }}} \right)_{p, {Y_i}}} = \frac{{RMT}}{{{{\left( {M - b\rho } \right)}^2}}} - a\frac{{2\rho {M^2} + \left( {{\delta _1} + {\delta _2}} \right)bM{\rho ^2}}}{{{{\left( {M + {\delta _1}b\rho } \right)}^2}{{\left( {M + {\delta _2}b\rho } \right)}^2}}} $ (6)
$ \begin{array}{*{20}{l}} {{{\left( {\frac{{\partial p}}{{\partial {\rho _i}}}} \right)}_{T, {\rho _{{\rm{j}} \ne {\rm{i}}}}}} = \frac{{MRT\left( {M + \left( {{b_i} - b} \right)\rho } \right)}}{{{M_i}{{\left( {M - b\rho } \right)}^2}}} - }\\ {\frac{{\rho M\left( {\sum\limits_{j = 1}^{ns} {x_j}{a_{ij}} + a} \right)}}{{{M_i}\left( {M + {\delta _1}b\rho } \right)\left( {M + {\delta _2}b\rho } \right)}} + }\\ {\frac{{a{b_i}{\rho ^2}M\left( {{\delta _1}M + 2{\delta _1}{\delta _2}\rho b + {\delta _2}M} \right)}}{{{M_i}{{\left( {M + {\delta _1}b\rho } \right)}^2}{{\left( {M + {\delta _2}b\rho } \right)}^2}}}} \end{array} $ (7)

式中M为摩尔质量,R为气体常数。当流体为混合物时,根据混合规则有

$ \frac{{\partial a}}{{\partial T}} = \sum\limits_{i = 1}^{ns} \sum\limits_{j = 1}^{ns} {X_i}{X_j}\frac{{\partial {a_{ij}}}}{{\partial T}} = \sum\limits_{i = 1}^{ns} \sum\limits_{j = 1}^{ns} {X_i}{X_j}\frac{{\partial \sqrt {{a_i}{a_j}} }}{{\partial T}} $ (8)

采用一般流体状态方程时,内能、焓、比热、声速等参数的求解不能再采用理想气体状态方程下的简单形式,需要根据热力学关系式进行重新推导。

$ e = {e_0} + \frac{1}{{\left( {{\delta _1} - {\delta _2}} \right)bM}}\left( {T\frac{{\partial a}}{{\partial T}} - a} \right){\rm{ln}}\left( {\frac{{M + {\delta _1}b\rho }}{{M + {\delta _2}b\rho }}} \right) $ (9)
$ h = {h_0} + \frac{1}{{\left( {{\delta _1} - {\delta _2}} \right)bM}}\left( {T\frac{{\partial a}}{{\partial T}} - a} \right){\rm{ln}}\left( {\frac{{M + {\delta _1}b\rho }}{{M + {\delta _2}b\rho }}} \right) $ (10)
$ {c_v} = {c_{v, 0}} + \frac{T}{{\left( {{\delta _1} - {\delta _2}} \right)bM}}\left( {\frac{{{\partial ^2}a}}{{\partial {T^2}}}} \right){\rm{ln}}\left( {\frac{{M + {\delta _1}b\rho }}{{M + {\delta _2}b\rho }}} \right) $ (11)
$ {c_p} = {c_v} + \frac{T}{{{\rho ^2}}}{\left( {\frac{{\partial p}}{{\partial T}}} \right)^2}_{{\rho _i}}/{\left( {\frac{{\partial p}}{{\partial \rho }}} \right)_{T, {Y_i}}} $ (12)
$ {c^2} = \frac{{{c_p}}}{{{c_v}}}\left[ {\frac{{RMT}}{{{{\left( {M - b\rho } \right)}^2}}} - a\frac{{2\rho {M^2} + \left( {{\delta _1} + {\delta _2}} \right)bM{\rho ^2}}}{{{{\left( {M + {\delta _1}b\rho } \right)}^2}{{\left( {M + {\delta _2}b\rho } \right)}^2}}}} \right] $ (13)

下标0代表理想气体在温度T下对应的值。

2.3 预处理方法

可压缩算法用于求解压缩性较强的流动问题时具有良好的收敛性,但对于低速问题,存在两点困难:一是对流项过小引起的压力项奇异问题,二是Jacobian矩阵的特征值差异过大导致的刚性问题。

在低速条件下,动量方程中的对流项明显小于压力项,压力的微小变化就会引起速度的明显改变,导致压力项奇异问题。可通过将压力p拆分成平均压力p0和脉动压力pg之和的方式解决[15]。由于对方程起作用的是压力梯度,拆分压力并不会导致方程形式的改变,却可以减小压力梯度的截断误差,从而消除压力项奇异问题。

$ \nabla p = \nabla \left( {{p_0} + {p_{\rm{g}}}} \right) = \nabla {p_{\rm{g}}} $ (14)

在低速问题中,流速和声速存在量级上的差距,导致控制方程的特征值差异过大,引起数值刚性。在定常求解过程中,NS方程中的时间项为虚拟项,不具有物理意义,其值最终会趋向于零,改变时间项的形式并不会影响结果。因此,采用“人工可压缩性方法”[16]对时间项进行预处理,通过引入预处理矩阵Γ,改变NS方程的特征值,消除低速问题下的数值刚性。

$ \frac{{\nabla \mathit{\boldsymbol{Q}}}}{{\nabla t}} = \frac{{\nabla \mathit{\boldsymbol{Q}}}}{{\nabla \mathit{\boldsymbol{U}}}}\frac{{\nabla \mathit{\boldsymbol{U}}}}{{\nabla t}} \to \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\frac{{\nabla U}}{{\nabla t}} $ (15)

Q为守恒变量{ρ, ρu, ρv, ρw, ρE, ρY1, …, ρYns-1},U为原始变量{p, u, v, w, h, Y1, …, Yns-1},预处理矩阵Γ采用Shuen[17]针对化学反应提出的形式。

$ \mathit{\boldsymbol{ \boldsymbol{\varGamma} }} = \left( {\begin{array}{*{20}{c}} \mathit{\Theta }&0&0&0&0&0& \cdots &0\\ {\mathit{\Theta }u}&\rho &0&0&0&0& \cdots &0\\ {\mathit{\Theta }v}&0&\rho &0&0&0& \cdots &0\\ {\mathit{\Theta }w}&0&0&\rho &0&0& \cdots &0\\ {\mathit{\Theta }H - 1}&{\rho u}&{\rho v}&{\rho w}&\rho &0& \cdots &0\\ {\mathit{\Theta }{Y_1}}&0&0&0&0&\rho&\cdots &0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots &{}& \vdots \\ {\mathit{\Theta }{Y_{ns - 1}}}&0&0&0&0&0& \cdots &\rho \end{array}} \right) $ (16)

Θ为换算系数,与∂ρ/∂p同量纲。

$ \mathit{\Theta } = \frac{1}{{U_{{\rm{ref}}}^2}} - \frac{1}{{{a^2}}} $ (17)

参考速度Uref=min[a, max(|V|, K|V|)]。K为系数,可根据不同问题取不同值,V取流场中的最大速度。引入预处理矩阵后,方程特征值由原来的{u, …, u, u±a}变为{u, …, u, u′±a′},其中

$ u' \pm a' = \frac{1}{2}\left[ {\left( {1 + {\varepsilon ^2}} \right)u \pm a\sqrt {{{\left( {1 - {\varepsilon ^2}} \right)}^2}{M^2} + 4{\varepsilon ^2}} } \right] $ (18)

a′为数值声速,ε为参考马赫数。当地马赫数大于1时,令ε=1,则u′±a′=u±a,此时的特征值与预处理前相同。当地马赫数很低时,ε→0,此时的特征值与速度u量级相当,从而消除刚性问题。

2.4 裂解反应模型

燃油裂解反应过程较复杂,侯凌云等[18]根据Ward提出的PPD(Proportional Product Distribution)模型[19],对燃油裂解反应机理进行了简化,采用C7H8来描述燃油裂解的液相平均产物。

$ \begin{array}{*{20}{l}} {{\rm{RP}} - 3 \to 0.16{{\rm{H}}_2} + {\rm{C}}{{\rm{H}}_4} + 0.58{{\rm{C}}_2}{{\rm{H}}_6} + 0.43{{\rm{C}}_2}{{\rm{H}}_4} + }\\ {0.42{{\rm{C}}_3}{{\rm{H}}_6} + 0.28{{\rm{C}}_3}{{\rm{H}}_8} + 0.25{{\rm{C}}_4}{{\rm{H}}_8} + 0.84{{\rm{C}}_7}{{\rm{H}}_8}} \end{array} $ (19)

PPD模型较为简单,可保证低裂解度下(小于20%)的模拟精度,由于没有考虑二次裂解的影响,高裂解度下的组分分布误差会增大。但二次裂解反应的吸热能力相对较弱[7],对温度分布的影响比一次裂解小得多,且本文主要关注可压缩算法在流动传热中的应用,对二次裂解产物成分不敏感,因此采用了相对简单的反应机理。反应速率由Arrhenius公式计算,指前因子取2.1 × 1015s-1,活化能取237.8kJmol-1

采用多重网格技术加速收敛时,由于化学反应速率与温度成非线性关系,由粗网格上物理量直接求化学反应源项可能产生较大的误差,导致粗细网格不匹配,出现数值不稳定。仅在最细网格上计算化学反应源项及雅克比矩阵,再插值至粗网格上使用。

2.5 物性计算方法

为考虑温度和压力对物性的影响,采用Chung方法[20]计算粘性系数。由于燃油裂解过程中存在较多的烃类化合物,由Brule-Starling方法[20]定义其中的系数值,以保证在较宽温度和压力范围内的精度。

$ \eta = {\eta ^{\rm{*}}}\frac{{36.344{{\left( {M{T_{\rm{c}}}} \right)}^{0.5}}}}{{V_{\rm{c}}^{2/3}}} $ (20)

式中Tc为临界温度,Vc为临界体积,M为分子量,η*值由临界参数及相应系数确定。

Chung引入压力的影响,提出了考虑高压影响的导热率计算方法。

$ \lambda = \frac{{31.2{\eta ^ \circ }\mathit{\Psi }}}{M}\left( {G_2^{ - 1} + {B_6}y} \right) + q{B_7}{y^2}T_{\rm{r}}^{1/2}{G_2} $ (21)

式中η°为低压下的粘性系数,Tr为对比温度,Ψqy为由临界参数确定的函数关系,BiGi为由偏心因子、对比偶极矩、缔合因子确定的函数关系。

针对低压下的双组元混合物,采用Chapman提出的扩散系数计算方法。并根据Dawson等提出的关系式,对扩散系数进行高压修正。

$ {D_{{\rm{AB}}}} = 1.858 \times {10^{ - 7}}\frac{{{T^{3/2}}\sqrt {1/{W_{\rm{A}}} + 1/{W_{\rm{B}}}} }}{{p\sigma _{{\rm{AB}}}^2{\mathit{\Omega }_{\rm{D}}}}} $ (22)
$ \begin{array}{l} {D_{{\rm{AB}}}}\rho = \\ {\left( {{D_{{\rm{AB}}}}\rho } \right)^\circ }\left( {1 + 0.053432{\rho _{\rm{r}}} - 0.030182\rho _{\rm{r}}^2 - 0.029725\rho _{\rm{r}}^3} \right) \end{array} $ (23)

式中σAB为分子间的碰撞直径,ΩD为碰撞积分,(DABρ)°代表低压下的扩散系数与密度的乘积,ρr为对比密度ρ/ρc

2.6 流热耦合边界

结构与流体之间的耦合传热采用分区计算、边界耦合的方式实现,即结构与流体各自求解自身的控制方程,仅在交界面处传递参数,交界面处需保证温度和热流密度的连续性。

$ {T_{{\rm{w}}, {\rm{solid}}}} = {T_{{\rm{w}}, {\rm{fluid}}}}\;\;\;{q_{{\rm{w}}, {\rm{solid}}}} = {q_{{\rm{w}}, {\rm{fluid}}}} $ (24)
3 计算结果分析 3.1 物性计算方法验证

RP-3燃油为一种包含近百种组分的混合物,通常采用几个主要组分进行替代,使其物性参数相近。目前存在多种燃油的替代模型,采用天津大学蒋榕培等[5]提出的四组分模型(25%正癸烷,50%正十二烷,12%正十三烷,13%正丁基环己烷)。通过本文采用的物性方法计算燃油在5MPa压力下的物性参数,并将其与文献结果进行对比,如图 1所示。可见,本文采用的物性计算方法与文献结果符合较好。

Fig. 1 Thermophysical properties variety of PR-3 with temperature(5MPa)
3.2 电加热管验证算例

天津大学蒋榕培等[5]通过开展电加热管试验,研究了燃油主动冷却过程。试验对镍基高温合金直圆管进行电加热,并由燃油进行冷却,试验条件如表 1所示。

Table 1 Operation conditions of electrically heated tube tests

由于试验涉及到的流动传热过程与超燃发动机的燃油冷却过程相近,对该试验进行数值模拟,以验证本文数值方法的有效性。图 2为管内流体的马赫数和密度的分布云图,图 3为管内流体的截面平均温度及燃油裂解率沿轴向的分布。

Fig. 2 Distributions of variables along the tube

Fig. 3 Variations of averaged variables along the tube

裂解率反映了燃油的裂解程度,其定义如式(25)所示,其中下标RP-3代表未反应的燃油质量,下标in代表进入计算域的燃油总质量。

$ \chi = \left( {1 - {m_{{\rm{RP}} - 3}}/{m_{{\rm{in}}}}} \right) \times 100{\rm{\% }} $ (25)

由图中结果可知,燃油裂解反应从0.4m处开始发生,0.5m后反应明显加快,对应燃油截面平均温度约800K。说明燃油在800K以下时,主要靠物理热沉吸热,燃油的温度基本呈线性增加。超过800K后,燃油开始快速裂解,化学热沉发挥主要作用,在较小的温度区间内由裂解反应吸收了大量热量,使燃油的温升速率显著降低。随着裂解率的增大,燃油的截面平均温度、裂解率的计算结果与试验结果偏差逐渐增大,金属管的出口处温度偏大约30K,裂解率偏大约7%。总体来说,计算结果与试验结果符合很好。

3.3 燃油冷却面板计算分析

电加热管构型较为简单,通过更接近真实情况的平板面板对燃油冷却过程进行研究。面板长350mm,宽90mm,厚6mm。由于燃油槽道并联容易产生流量分配不均的问题,因此采用蛇形串联布置,槽道总长超过5m,截面为矩形,宽2mm,高1mm。研究了两种不同的进出油方式,分别为一侧进油一侧出油和两侧进油中间出油,如图 4所示。面板一侧加载1MW/m2的热流,加热范围350mm×75mm。

Fig. 4 Cooling channel arrangement

针对两侧进油方式,计算了不同燃油流量下的冷却面板温度分布,燃油入口压力5MPa。如图 5所示,冷却面板温度分布呈现两侧低中间高的趋势。随着燃油流量的增大,面板的温度明显下降。四种不同流量下,燃油出口最高温度分别达到了951K,884K,787K和694K。对于RP-3燃油,在温度超过900K以后,会发生明显的结焦现象,严重时可能导致冷却通道堵塞,燃油流量减少,面板温度迅速升高。由于计算过程中没有采用结焦模型,无法反映结焦对传热过程的影响。虽然计算给出的面板最高温度均没有超过1100K,仍在高温合金的材料许可范围内,但实际上燃油流量为7.3g/s时,燃油最高已经达到951K,冷却通道内将会发生严重的结焦,面板可能局部温度过高而导致防热失效。因此,在当前的计算条件下,燃油流量大于10g/s时,可以保证防热安全。

Fig. 5 Distribution of temperature under different mass flow rates

计算了两种不同进出油方式下的冷却面板温度分布,总流量10g/s,两侧进油方式每个槽道各5g/s,如图 6所示。面板的温度分布与燃油的温度分布保持一致,入口处温度低,出口处温度高。两种进出油方式的面板和燃油最高温度都在许用的范围内,但单侧进油情况下的面板最大温差约380K,比两侧进油情况高了近150K。显然,两侧进油情况下的面板均温性更好,热应力更低。

Fig. 6 Distribution of temperature under different cooling channel arrangements
4 结论

本文通过开展数值研究,获得如下结论:

(1)采用的物理模型及数值方法在燃油冷却过程的数值模拟中是准确可靠的。

(2)800K为RP-3燃油的重要温度节点,低于800K时主要靠物理热沉吸热,高于800K时主要靠裂解产生的化学热沉吸热。

(3)在1MW/m2热流下,尺寸350mm×90mm×6mm的平板面板需要至少10g/s的RP-3燃油才能保证防热安全。

(4)对于串联式槽道布置,在燃油流量相同情况下,两侧进油方式比单侧进油方式最大温差低约150K,均温性更好,更有利于降低结构热应力。

在今后的研究工作中,还应进一步分析入口流态、通道流阻等因素的影响,以给出满足实际工程需求的冷却面板设计。

参考文献
[1]
范学军, 俞刚. 大庆RP-3航空煤油热物性分析[J]. 推进技术, 2006, 27(2): 187-192. (FAN Xue-jun, YU Gang. Analysis of Thermophysical Properties of Daqing RP-3 Aviation Kerosene[J]. Journal of Propulsion Technology, 2006, 27(2): 187-192. DOI:10.3321/j.issn:1001-4055.2006.02.021) (0)
[2]
王夕.超临界压力吸热型碳氢燃料热裂解及传热特性研究[D].北京: 清华大学, 2013. http://cdmd.cnki.com.cn/article/cdmd-10003-1014021223.htm (0)
[3]
赵国柱, 宋文艳, 张若凌, 等. 超临界压力下正十烷流动传热的数值模拟[J]. 推进技术, 2014, 35(4): 537-543. (ZHAO Guo-zhu, SONG Wen-yan, ZHANG Ruo-ling, et al. Numerical Simulation on Flow and Heat Transfer of n-Decane under Supercritical Pressure[J]. Journal of Propulsion Technology, 2014, 35(4): 537-543.) (0)
[4]
赵国柱, 宋文艳, 张若凌. 超临界压力下RP-3航空煤油吸热裂解反应的数值研究[J]. 航空学报, 2014, 35(6): 1513-1521. (0)
[5]
Jiang R P, Liu G Z, Wang X Q, et al. Thermal Cracking of Hydrocarbon Aviation Fuels in Regenerative Cooling Microchannels[J]. Energy and Fuels, 2013, 27(5): 2563-2577. DOI:10.1021/ef400367n (0)
[6]
Zhang R L, Zhao G Z, Le J L, et al. Numerical Study on Heat Transfer of Hydrocarbon Fuel with Thermal Cracking[R]. AIAA 2015-3621. (0)
[7]
张强强.微细通道内碳氢燃料传热与裂解过程的实验研究与CFD模拟[D].天津: 天津大学, 2014. http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=D636808 (0)
[8]
Bao W, Zhang S L, Qin J, et al. Numerical Analysis of Flowing Cracked Hydrocarbon Fuel inside Cooling Channels in View of Thermal Management[J]. Energy, 2014, 67(4): 149-161. (0)
[9]
Zhang S L, Qin J, Xie K L, et al. Thermal Behavior inside Scramjet Cooling Channels at Different Channel Aspect Ratios[J]. Journal of Propulsion and Power, 2016, 32(1): 57-70. DOI:10.2514/1.B35563 (0)
[10]
Menter F R. Two Equation Eddy Viscosity Turbulence Models for Engineering Applications[J]. AIAA Journal, 2012, 32(8): 1598-1605. (0)
[11]
Yoon S, Jameson A. Lower-Upper Symmetric Gauss-Seidel Method for the Euler and Navier-Stokes Equations[J]. AIAA Journal, 1988, 26(9): 1025-1026. DOI:10.2514/3.10007 (0)
[12]
Edwards J R, Franklin R K, Liou M S. Low-Diffusion Flux-Splitting Methods for Real Fluid Flows at All Speeds[R]. AIAA 99-3327. (0)
[13]
Kim S K, Choi H S, Kim Y. Thermodynamic Modeling Based on a Generalized Cubic Equation of State for Kerosene/LOX Rocket Combustion[J]. Combustion and Flame, 2012, 159(3): 1351-1365. DOI:10.1016/j.combustflame.2011.10.008 (0)
[14]
Cismondi M, Mollerup J. Development and Application of a Three-Parameter RK-PR Equation of State[J]. Fluid Phase Equilibria, 2005, 232(1): 74-89. (0)
[15]
Meng H, Yang V. A Unified Treatment of General Fluid Thermodynamics and its Application to a Preconditioning Scheme[J]. Journal of Computational Physics, 2003, 189(1): 277-304. DOI:10.1016/S0021-9991(03)00211-0 (0)
[16]
Chorin A J. A Numerical Method for Solving Incompressible Viscous Flow Problem[J]. Journal of Computational Physics, 1967, 2(1): 12-26. DOI:10.1016/0021-9991(67)90037-X (0)
[17]
Shuen J S, Chen K H, Choi Y. A Coupled Implicit Method for Chemical Non-Equilibrium Flows at All Speeds[J]. Journal of Computational Physics, 1993, 106(2): 306-318. DOI:10.1016/S0021-9991(83)71110-1 (0)
[18]
侯凌云, 董宁, 孙大鹏. 矩形冷却槽道内煤油热裂解过程数值研究[J]. 推进技术, 2014, 35(1): 128-132. (HOU Ling-yun, DONG Ning, SUN Da-peng. Numerical Study on Thermal Cracking Process of Kerosene in a Rectangular Cooling Channel[J]. Journal of Propulsion Technology, 2014, 35(1): 128-132.) (0)
[19]
Ward T A, Ervin J S, Striebich R C, et al. Simulations of Flowing Mildly-Cracked Normal Alkanes Incorporating Proportional Product Distributions[J]. Journal of Propulsion and Power, 2004, 20(3): 394-402. DOI:10.2514/1.10380 (0)
[20]
童景山. 流体热物性学基本理论与计算[M]. 北京: 中国石化出版社, 2008. (0)