基于热解动力学炭/酚醛燃气舵流热耦合数值研究

薛海峰，陈 雄，郑 健，周长省

(南京理工大学 机械工程学院，南京 210094)



基于热解动力学炭/酚醛燃气舵流热耦合数值研究

薛海峰，陈 雄，郑 健，周长省

(南京理工大学 机械工程学院，南京 210094)

针对炭/酚醛燃气舵体积烧蚀问题，在Fluent平台上利用UDF二次开发进行了二维非定常流热耦合数值研究。对几何建模、材料变热物性模型及边界条件等问题进行了详尽的描述，并选取了合适的计算模型。对不同舵偏角下燃气舵温度分布、材料密度及边界热流密度等参数进行了分析研究。计算结果表明，燃气舵前缘一直是体积烧蚀最严重区域，随着舵偏角的增大，迎风面体积烧蚀越为严重；由于炭/酚醛材料的特殊性，随着工作时间的推进，从边界进入燃气舵内部热流密度逐渐降低，趋于一个稳定值。研究方法及结论可用于炭化烧蚀类复合材料燃气舵热分析研究。

炭/酚醛；燃气舵；流热耦合；体积烧蚀

0 引言

以燃气舵为控制方式的推力矢量控制系统，目前被广泛用于各型导弹中[1]。燃气舵在整个工作过程中，都处于高温超声速燃气射流氛围内，其工作环境极其恶劣。为了克服燃气射流对燃气的烧蚀影响，当前燃气舵使用的主流材料为钨渗铜，但钨渗铜材料密度较大，降低了弹箭的工作效率[2-3]。

目前,国内外学者已经对炭纤维增强复合材料燃气舵进行了相关实验研究。Kumar Suresh等[4]对碳布采用穿刺工艺制成类三维编织增强体的碳/碳化硅燃气舵，并以含铝颗粒的高能复合装药固体火箭发动机对该燃气舵进行地面烧蚀试验；Chen Bo等[2]对针刺C/C复合材料机加工制成的楔形块进行了固体火箭发动机燃气射流烧蚀试验及相关流动仿真研究；Bansard S，Plouvier S等[5-6]对炭/酚醛复合材料进行了含液态铝高温气流下的烧蚀实验研究。相关试验研究表明，炭纤维增强复合材料具有很好的绝热和耐烧蚀性能。由于燃气舵的工作环境为高温超声速射流，且含有大量烟气，对燃气舵表面温度及内部温度场的实验方法获取存在一定难度，借助有效的数值仿真手段，可得到燃气舵在工作过程中表面及内部温度场分布。

对于燃气舵传热过程的研究，绝大部分学者[7-10]都是通过单向耦合方式来处理。具体而言，即通过绕流流场的定常求解来获取固体区域的边界条件，并将该边界条件用于燃气舵温度场的非定常求解。该方法避免了绕流流场的非定常求解，提高了解决问题的效率。事实上，燃气舵内部温度场的分布同样会对绕流流场产生一定的影响，尤其对炭/酚醛材料而言，由于材料物性参数不断在发生变化，且热解吸热过程通过壁面也影响着流场结构。

本文通过紧耦合的方式，对不同舵偏角的二维燃气舵表面及内部温度场分布进行了仿真研究。通过UDF编写发动机燃气及材料物性随工作时间及温度的改变过程，在燃气舵内部网格添加能量源项，来模拟酚醛热解的吸热过程。

1 物理模型和计算方法

1.1 基本假设

本文研究基于以下假设：

(1)火箭发动机射流处于化学平衡状态，且为理想气体。

(2)酚醛树脂热解产物组分为甲烷(CH4)、乙烯(C2H4)、乙炔(C2H2)和苯蒸气(C6H6)[11]，且各组分摩尔分数为常数。

(3)不考虑热解气体从燃气舵内部逸出过程，由于热解气体逸出导致的热阻塞效应以能量源项的形式给出；忽略燃气逸出在壁面的引射效应。

(4)研究重点在于对燃气舵内部温度场以及体积烧蚀[12]，忽略由于热解导致的体积变化及热膨胀；不考虑热化学烧蚀以及机械剥蚀，因而不存在燃气舵表面退移现象。

(5)忽略辐射换热，忽略重力等体积力的影响。

1.2 计算模型及工况

本文在Fluent平台上，使用UDF修正了火箭发动机燃气热物性参数随温度的变化，并在燃气舵内部网格上，定义了能量源项来模拟炭/酚醛热解吸热过程及热解气体逸出所携带的能量。考虑到在超声速射流中燃气舵的存在，使得全场流速变化较大，本文使用k-ωSST两方程湍流模型。

1.2.1 炭/酚醛热解炭化模型

炭/酚醛复合材料的耐高温机制已经研究的非常明确，按照烧蚀机理将其划分为热解炭化类材料。根据炭/酚醛复合材料的传热烧蚀机理，通用的3层模型如图1所示，包括Ⅰ炭化层、Ⅱ热解层和Ⅲ原始材料层。罗永康等[13]的研究表明，在2 800 K温度以下炭化形成的酚醛树脂碳主要成分为无定形碳(Free Carbon)。本文假设炭化层中的主要构成为无定形碳、炭纤维以及酚醛树脂热解形成的孔隙；热解层中的主要构成为酚醛树脂、部分酚醛树脂热解形成的无定形碳、炭纤维及热解形成的孔隙；原始材料层中主要成分为酚醛树脂和炭纤维。

图1 炭/酚醛复合材料3层模型示意图Fig.1 Schematic of three layer model about carbon-phenolic

1.2.2 物理模型及边界条件

燃气舵绕流流场是一个三维流动，由于燃气舵内部传热是一个非定常过程，对计算机的硬件要求较高。文献[14-15]指出，在满足仿真要求的前提下，将三维模型简化为二维模型是可行的。本文对舵偏角为0°、5°、10°、15°和20°等5种情况下燃气舵的流热耦合进行了数值仿真研究。以10°舵偏角为例，图2给出了仿真模型的结构示意图以及边界条件；为了方便分析说明，图2右下角给出了燃气舵放大图及舵面参考系说明。边界条件参数设置为：(1)压力远场，静压101 325 Pa，静温300 K；(2)压力入口，根据发动机内弹道参数给出总压8 MPa，静压120 715 Pa，总温2 766 K；(3)压力出口，反压101 325 Pa，静温300 K；(4)炭/酚醛燃气舵作为固体区域参与能量传递计算，舵面前缘、后缘、迎风面和背风面均设为耦合壁面。整个计算模型均划分为结构网格，固体区域网格单元数为9 346，流体区域网格单元数为41 270。

图2 仿真模型结构示意图Fig.2 Schematic of simulation model

1.2.3 物性参数及源项处理

(1)炭/酚醛密度模型

基于“多组分模型[16]”，给出炭/酚醛原始材料密度表达式：

ρv=ΓρA+(1-Γ)ρB

(1)

在热解过程中，以阿累尼乌斯方程来表征酚醛树脂的热解速率[11]，材料密度的变化过程为

(2)

在当前时间步内积分式(2)，且在该时间步内视温度T为常量，有

(3)

文中炭化层密度ρch取决于材料初始状态的配方及酚醛树脂的成碳率，由式(4)给出：

ρch=εcharΓρA+(1-Γ)ρB

(4)

(2)炭/酚醛比热容模型

由于在Fluent平台上，材料的热容只能定义为与温度相关的函数。因而在炭/酚醛比定压热容的修正上，根据文献[17]采用关于温度的多项式拟合如下：

cp=a+bT+cT2+dT3+eT4+fT5+gT6

(5)

当300≤T<1 922时，a=-1 308.33，b= 11.77，c= -1.72×10-2，d= 9.86×10-6，e= -7.15×10-11，f= -1.83×10-12，g= 4.61×10-16；当1 922≤T<5 000时，a= 2 023.85，b= 5.52×10-2，c～f均为0。

(3)炭/酚醛热导率模型

由于炭化层和热解层中存在大量的孔隙，基于多孔介质传热理论[18]中“有效导热系数法”，给出炭/酚醛在不同状态下的热导率表达方程：

k=ΓAkA+ΓBkB+ΓCkC+ΦkG

(6)

(7)

(4)热解气体和发动机燃气热物性模型

炭/酚醛燃气舵在工作过程中，内部温度变化较大，这会导致热解气体导热系数及比空压热容变化较大，为保证数值计算的准确性，需要对其进行估算。表1给出了酚醛树脂高温热解产物分布[11]。

表1 热解产物摩尔分数Table1 Mole fraction of pyrolysis products

忽略压强变化对热解气体热物性的影响。查阅JANAF表，可获得表1中各组分在不同温度下比定压热容(J/(kg·K))，进行分段多项式拟合如下形式：

cp=a+bT+cT2+dT3+eT4

(8)

具体结果见表2。

表2 各组分比定压热容多段拟合系数Table2 Coefficient of polynomial fitting about specific heat of each component

根据多原子气体的Eucken A关系式及表2中比空压热容函数，可得出各热解组分热导率以及粘性系数随温度变化关系。

根据最小自由能法计算出发动机燃气组分分布，具体见表3，发动机燃气的等效分子量为24.627。由于本文研究不考虑发动机燃气与舵面的热化学反应，为了提高计算效率，同样参照热解气体热物性方法获取燃气的比空压热容、热导率及粘性系数。其中，各组分比空压热容多项式见表2。

表3 燃气组分摩尔分数Table3 Mole fraction of rocket gas

(5)能量源项

能量源项S(W/m3)包含了酚醛树脂热解潜热及热解气体逸出所携带能量的等效部分：

S=Spy+Sen

(9)

其中

式中Spy为酚醛树脂热解能量源项；Sen为热解气体逸出携带的能量源项；hA为酚醛树脂热解潜热；hg为热解气体显焓。

上述对炭/酚醛热物性参数以及源项处理方式的具体过程可参考文献[19]，且通过该文献中实验验证了模型的准确性。

2 计算结果及分析

使用上述方法及参数，对炭/酚醛燃气舵在不同偏角姿态下进行非定常流热耦合数值模拟。炭/酚醛相关参数如下：酚醛树脂密度1 186 kg/m3，炭纤维密度为1 800 kg/m3，酚醛树脂的体积分数为0.6，酚醛树脂的成碳率为0.5，酚醛树脂完全炭化后的真实密度为1 500 kg/m3；热解反应临界温度为573 K，反应指前因子为185 000 s-1，活化能为100 810 J/mol[11]，热解潜热为420 000 J/kg[20]；燃气舵内部初始温度与环境温度保持一致为300 K。

2.1 不同舵偏角的影响

图3为t=2.0 s时刻，燃气舵在舵偏角为0°、10°、20°等3种情况下，舵体内部及其附近燃气绕流温度云图。在超声速燃气射流环境下，舵面前缘形成了弓形激波，通过激波的燃气温度骤增，燃气舵前缘表面温度基本为燃气总温。而在舵面后缘靠近壁面附近，由于该处气流分离形成驻涡，流动基本滞止，同样形成一个局部高温区。一旦燃气舵偏转，在迎风面一侧，斜激波波角要比背风面一侧激波波角大，因而迎风面一侧波后流场温度要高于背风面一侧温度。

(a)α=0° (b)α=10° (c)α=20°

图4和图5分别为t=2.0 s时，在不同舵偏角下燃气舵表面温度分布曲线和A-A面(即x=25 mm处)燃气舵内部y方向温度分布。由于流动的滞止作用，无论舵面偏转角度如何，燃气舵前缘壁面温度基本接近于燃气总温，为2 583 K。当流动越过了前缘与舵体侧面的圆弧过渡段后，舵面温度迅速下降。从该点往后，不同舵偏角导致的迎风面和背风面壁面温度差异较高。α=0°时，流动对称，因而侧面温度保持一致；侧壁面温度随着x的增大而降低。随着舵偏角的增大，迎风面壁面温度逐渐升高。这主要是由于该侧面流场温度较高，且在迎风面处来流方向与侧壁面存在一定的夹角，流动运动到此处发生滞止所引发的。反之，随着舵偏角的增加，背风面壁面温度逐渐降低，该处舵面并未受到燃气的直接冲击。当α=20°时，燃气舵位于迎风面侧面温度最高约为2 300 K，背风面最低温度位于燃气舵侧面末端，仅有1 500 K。

图5中虚线表示炭/酚醛热解临界温度573 K，超过该温度材料发生热解。当α=20°时，迎风面燃气舵开始发生热解区域距离舵面1.350 4 mm，背风面热解区域距离舵面0.906 1 mm。

图4 t=2.0 s时燃气舵表面温度Fig.4 Surface temperature of jet-vane when t=2.0 s

图5 t=2.0 s时，x=25 mm处燃气舵内部温度Fig.5 Internal temperature of jet-vane at x=25 mm when t=2.0 s

由于燃气舵在工作过程中一直处于高温超声速燃气射流氛围中，舵面炭化层厚度一旦超越某个临界值，在燃气射流的冲刷下会发生机械剥蚀[21]，燃气舵由于表面剥蚀导致厚度变薄，会存在断裂失效的可能。在本研究的数学模型中，密度作为燃气舵内部分层模型的判据，表征了燃气舵体积烧蚀程度。图6和图7分别为t=2.0 s时，不同舵偏角下燃气舵内部密度分布以及x=25 mm处炭化层厚度变化曲线。由图7可见，无论舵面如何偏转，燃气舵前缘的炭化层厚度最高；在燃气舵后缘，同样由于局部高温区的作用，也存在一个局部高炭化层厚度区域。当α=0°时，燃气舵两侧炭化层厚度保持一致，为0.454 0 mm；随着舵偏角的增加，燃气舵迎风面炭化层厚度逐渐增加，其中α=5°和10°时，炭化层厚度为0.564 9 mm，当α=15°和20°时，炭化层厚度增大到0.697 3 mm；而背风面炭化层厚度降低几乎以线性规律从0.454 0 mm降低到0.139 6 mm。

2.2 燃气舵传热烧蚀研究

以舵偏角为20°为例，分析了不同时刻点相关数据。图8为不同时刻点燃气舵在20°舵偏角状态下舵面温度分布曲线。在发动机工作初始时刻，燃气舵表面温度迅速升高。当t=0.01 s时，燃气舵前缘壁面温度已经达到1 394 K，在迎风面壁面温度也达到了900 K，此刻燃气舵前缘和迎风面均已开始发生热解；而背风面壁面温度仅不到500 K，并未发生热解。随着时间的推移，燃气舵壁面局部温度上升速度逐渐降低。这主要是由于酚醛树脂的热解吸热、热解气体逸出过程中携带了部分热量，以及材料本身的热容所造成的。

图7 t=2 s时，x=25 mm处炭化层厚度Fig.7 Thickness of charring layer at x=25 mm when t=2.0 s

图9为不同时刻点燃气舵表面热流密度分布曲线。当t=0.01s时，前缘流动滞止区域热流密度高达15.9 MW/m2，在迎风面和背风面热流密度分别达到了7 MW/m2和2 MW/m2左右；这主要是由于燃气舵工作初始时刻温度仅为300 K，高温燃气冲击到较低温度的舵面会产生很高的热流密度。当燃气舵工作到一定时间之后，随着表面温度的升高，热流密度逐渐下降。图10为不同时刻点x=25mm处燃气舵侧面炭化层厚度变化曲线。在t=0.5 s时，舵偏角为10°和20°燃气舵背风壁面处并未出现炭化现象，而此时迎风壁面处炭化层厚度已经分别达到了0.084 4 mm和0.1396 mm；而在20°舵偏角工况下，t=1.0 s时背风壁面依然没有出现炭化现象，这说明舵偏角对炭化层的生成有很大的影响。随着时间的推移，炭化层厚度逐渐增大，而迎风壁面炭化层厚度增长速率要高于背风面。

图8 α=20°时燃气舵表面温度分布Fig.8 Surface temperature of jet-vane at α=20°

(a)t=0.01 s

(b)t=0.5 s,1.0 s,1.5 s,2.0 s

2.3 体积烧蚀对壁面温度的影响

作为对比，对相同来流条件且不考虑体积烧蚀时，0°、10°和20°偏角下的燃气舵进行了流热耦合计算。并对比了2种耦合方式下燃气舵壁面以及内部温度分布。

图11为不考虑体积烧蚀前提下，t=2.0 s时燃气舵内部及绕流温度场。与图3对比，很明显由于没有考虑体积烧蚀的影响，燃气舵内部能量扩散较快，壁面内部高温区厚度近乎为图3中的2倍。

为了更详细地说明不考虑体积烧蚀条件下燃气舵壁面以及内部温度分布，图12给出了不同时刻点两者壁面温度对比曲线。

图10 α=20°，x=25 mm处炭化层厚度Fig.10 Thickness of charring layer at x=25 mm while α=20°

从图12(a)可发现，当t=0.01 s时，由于不考虑体积烧蚀，壁面温度要比体积烧蚀模型温度高约200 K；随着时间的推进，2个不同模型导致的壁面温度差距越来越大。这意味着由于热解烧蚀的存在，有效阻止了能量向燃气舵内部的传递。当舵偏角变大时，不考虑体积烧蚀的情况下，燃气舵迎风面和背风面温差逐渐增大。

4 结论

(1)炭/酚醛材料由于其较低的热导率，且酚醛树脂在热解过程中存在吸热现象，有效阻碍了热量向燃气舵内部的传递。

(2)在工作起始阶段，炭/酚醛燃气舵前缘温度迅速上升，迎风面温度升幅略低，背风面温度升幅最小；随着工作时间的推进，舵面当地温度上升速率逐渐降低。燃气舵前缘流动滞止区域始终是全舵面温度最高点，且表面温度沿着燃气流动方向逐渐降低，在燃气舵后缘存在一个局部二次高温区域。一旦燃气舵发生偏转，迎风面温度要高于背风面；随着舵偏角的增加，迎风面和背风面温差逐渐变大。

(3)燃气舵前缘是材料体积烧蚀最严重的区域，在后缘存在一个局部体积烧蚀严重区域。舵偏角对燃气舵侧面体积烧蚀有很大的影响，舵面发生偏转之后，迎风面体积烧蚀要比背风面体积烧蚀严重。

(4)本文给出的结果是在不考虑燃气舵表面退移的假设前提下得到的计算值。而在实际的研究应该是包含由于热化学烧蚀和机械剥蚀导致的炭/酚醛燃气舵壁面退移的流固热耦合问题，这是今后研究中的重点方向。

(a)α=0° (b)α=10° (c)α=20°

(a)α=0° (b)α=10° (c)α=20°

[1] 刘志珩.固体火箭燃气舵气动设计研究[J].导弹与航天运载技术,1995 (4):9-17.

[2] Chen B,Zhang L T,Cheng L F,et al.Erosion Resistance of Needled Carbon/Carbon Composites Exposed to Solid Rocket Motor Plumes[J].Carbon,2009,47(6):1474-1479.

[3] 刘丽丽,李克智,李贺军.基于有限元的C/C燃气舵振动特性[J].玻璃钢/复合材料,2011,216(1):12-15.

[4] Kumar S,Kumar A,Sampath K,et al.Fabrication and Erosion Studies of C-SiC Composite Jet Vanes in Solid Rocket Motor Exhaust[J].Journal of the European Ceramic Society,2011,31(13):2425-2431.

[5] Bansard S,Plouvier S,et al.Experimental Simulation of Thermo-Mechanical Ablation of Carbon/Phenolic Composite Under the Impact of Liquid Alumina Particles[J].High Temperature Material Processes,2005,9(3):431-441.

[6] Bansard S,Legros E,Vardelle M,et al.Diagnostics of Two plasma jet flows close to a carbon phenolic composite target using molecular emission spectroscopy[J].High Temperature Material Processes,2005,9(3):415-429.

[7] Yu M S,Lee J W,Cho H H,et al.Numerical Study on A Thermal Response of the Jet Vane System in A Rocket Nozzle[C]//42nd AIAA Conference.Reno,NV.2004.

[8] 董晓芳.固体火箭发动机燃气舵热分析研究[D].西北工业大学,2005.

[9] 李军,常见虎,周长省,等.推力矢量燃气舵三维气-固两相流的数值分析[J].南京理工大学学报:自然科学版,2008,32(5):565-569.

[10] Yu M S,Cho H H,Hwang K Y,et al.A Study on A Surface Ablation of the Jet Vane System in A Rocket Nozzle[C]//37th AIAA Conference.2004:28.

[11] 马伟.酚醛树脂的热解研究[D].重庆：重庆大学,2007.

[12] 易法军,梁军.防热复合材料的烧蚀机理与模型研究[J].固体火箭技术,2000,23(4):48-56.

[13] 罗永康,彭维舟,王为民.烧蚀复合材料用酚醛树脂研究[J].宇航材料工艺,1988(5):7.

[14] Cho H H,Kim B G,et al.Analysis of Particle Laden Flows around A Jet Vane in A solid Rocket motor[C]//37th AIAA Conference.July,2001.

[15] 曹熙炜,刘宇,谢侃,等.一种特型燃气舵数值模拟分析[J].固体火箭技术,2011,34(1):5-8.

[16] Chen Y K,Milos F S.Two-Dimensional Implicit Thermal Response and Ablation Program for Charring Materials[J].Journal of Spacecraft and Rockets,2001,38(4):473-481.

[17] Potts R L.Application of Integral Methods to Ablation Charring Erosion,A Review[J].Journal of Spacecraft and Rockets,1995,32(2):200-209.

[18] 林瑞泰.多孔介质传热传质引论[M].科学出版社,1995.

[19] 薛海峰,陈雄,周长省.基于热解过程的变热物性炭/酚醛能量扩散数值研究[J].固体火箭技术,38(1):130-135.

[20] 徐晓亮.热防护机理与烧蚀钝体绕流的涡方法研究[D].北京：北京交通大学,2011.

[21] DiCristina V,Howey D,et al.Thermomechanical Erosion of Ablative Plastic Composites[R].Avco Systems Div Wilmington Mass,1971.

(编辑：吕耀辉)

Numerical research on flow-thermal coupling of carbon-phenolic jet-vane based on pyrolysis kinetics

XUE Hai-feng,CHEN Xiong,ZHENG Jian,ZHOU Chang-sheng

(School of Mechanical Engineering,Nanjing University of Science and Technology,Nanjing 210094,China)

To solve the problem of volumetric ablation about carbon-phenolic jet-vane, a numerical research on two-dimensional unsteady flow-thermal coupling was done using secondary development on fluent platform.Geometric model,variable thermal properties of carbon-phenolic as well as boundary conditions were described carefully,and a suitable numerical model was selected.Temperature distribution of jet-vane,material density and surface heat flux of jet-vane under different deflection angles of control surface were analyzed.The simulated results show that leading edge is always the most serious volumetric ablative area of jet-vane.As the deflection angle of jet-vane increases,volumetric ablation of upwind side becomes much more serious.Due to the particularity of carbon-phenolic,along with the working time, heat flux into jet-vane from boundary is reduced gradually,and tends to be a stable value.Research method and conclusions can be used in study of jet-vane made of charring materials.

carbon-phenolic;jet vane;flow-thermal coupling;volumetric ablation

2014-09-27；

：2014-11-17。

薛海峰(1986—)，男，博士生，研究领域为固体火箭发动机热防护。E-mail：liangwangongli@163.com

V421.6+2

A

1006-2793(2015)04-0503-07

10.7673/j.issn.1006-2793.2015.04.010