井下火箭发射燃气对消音层冲蚀损伤模型及特性

李良，冯永保，马长林，夏文龙

(火箭军工程大学 导弹工程学院, 陕西 西安 710025)

0 引言

地下井火箭发射时，发动机喷出的高温、高速燃气流与周围相对静止空气介质快速混合形成强烈的湍流脉动，产生巨大噪音，不仅会对井内设施造成危害，也会降低箭载仪器的可靠性[1-2]。为降低井内噪音，消除噪音给火箭发射带来的安全隐患，在地下井内壁周围建有消音层，常见的材料有薄层微孔钢板和玻璃丝棉。火箭发动机点火后，由于地下井内部空间有限，且排烟道狭长，火箭出井过程中井内环境恶劣，存在回流引射效应与压力脉冲等问题，这些因素的综合作用会导致消音层发生烧蚀损伤，甚至影响发射安全[3-4]。因此，研究火箭发射出井过程中燃气流对消音层的冲蚀损伤规律，对消音层的后期修复以及地下井的结构优化具有重要意义。

由于问题的特殊性，对地下井消音层燃气烧蚀问题的研究国内外参考文献较少，但其根本属于燃气流场中的热-流体-固体多场耦合问题，涉及多场耦合方法、耦合界面数据传递、材料烧蚀动边界处理等技术，国内外学者都进行了大量研究。Bendiksen[5]基于任意拉格朗日-欧拉(ALE)方法，通过建立全耦合模型，对气弹问题进行了数值仿真计算；Bathe等[6]通过引入有限体积法，将全耦合数值仿真方法扩展到求解高速流场中的多场耦合问题；吴永海等[7]针对速射火炮身管的温度场计算，建立了统一的流体与固体(简称流固)传热数学模型，并用计算流体力学(CFD)方法得到了身管系统的瞬时温度历程；Degroote等[8]在分区拟牛顿迭代法基础上，提出了一种新的紧耦合分区迭代方法，提高了计算效率；张艳岗等[9]针对柴油机气缸盖工作数值仿真模拟，提出了一种基于整场离散的热-流体-固体多物理场耦合分析方法，得到了气缸盖的温度场；汪超台等[10]针对发动机排气涡轮在高温及气动载荷、离心载荷下的数值仿真模拟，联合流体动力学和有限元分析进行顺序流固耦合传热分析，计算出排气涡轮的温度场，并进行了热应力分析；金文奇[11]建立了火炮全膛烧蚀磨损内弹道模型，导出模型解算方法，仿真表明该模型和诸元解算方法是正确的，可用于烧蚀磨损火炮身管的寿命预测。

由于地下井内燃气流场极为复杂，在对消音层燃气冲蚀损伤数值仿真模拟时，若采用双向热-流体-固体耦合方法，同时对燃气流场、消音层温度场与变形场进行计算，会导致计算量庞大、不易收敛等问题。考虑到消音层冲蚀损伤对燃气流场的影响十分有限，本文忽略消音层烧蚀对燃气流场的影响，采用单向热-流体-固体耦合方法，利用前期得到的地下井热发射燃气流场数据，通过编写程序实现耦合交界面上的数据传递，在有限元分析软件Abaqus中建立有限元分析模型，分析燃气流对消音层的冲蚀损伤情况，并分析地下井直径、燃气射流初始冲击高度及发动机建压速率3个因素对玻璃丝棉烧蚀损伤的影响。

1 玻璃丝棉消音层冲蚀损伤数学模型

火箭出井过程中，井内传热过程极为复杂，包括燃气流不同温度区域的热传导，耦合交界面上燃气流与消音层内壁面之间的强迫对流换热和辐射换热、消音层内部的热传导等过程。考虑研究问题的复杂性，为简化模型，消音层传热计算时，假设燃气和消音层壁面之间只存在对流换热，暂时忽略辐射传热的影响。

1.1 玻璃丝棉相变传热数学模型

靠近燃气流一侧的玻璃纤维在高温燃气作用下不断升温，达到熔点时开始熔化，形成以熔融SiO2为主要成分的高黏性液态层，表层以下玻璃丝棉温升远不如表层，其受热情况符合半无限大物体的传热特点。玻璃丝棉相变传热属于具有移动边界的传热过程，以内壁面某一点为坐标原点，沿厚度方向x轴建立一维坐标系，相界面的位置随时间推移不断向固相区移动，设为s(t)，相界面处有相同的温度即相变温度，相界面为固相和液相区域共同的移动边界。根据导热基本定律与能量守恒定律可知：

玻璃丝棉固相区导热微分方程

(1)

初始条件

Ts(x,t)|t=0=Ti，

(2)

边界条件

Ts(x,t)|x=∞=Ti；

(3)

玻璃丝棉液相区导热微分方程

(4)

初始条件

Tl(x,t)|x=0=Tw，

(5)

边界条件

Ts(x,t)|x=s(t)=Tl(x,t)|x=s(t)=Tp，

(6)

(7)

式中：Ts(x,t)、Tl(x,t)分别为固相和液相温度；as、al分别为固相和液相导温系数；Ti为固相初始温度；Tw为壁面温度；Tp为相变温度；λs为固相热传导系数；λl为液相热传导系数；ρs为固相密度；L为熔化潜热。

1.2 玻璃丝棉烧蚀数学模型

玻璃丝棉烧蚀的特点是，随着烧蚀过程的发生，玻璃丝棉逐渐变薄，熔化形成的液态层在燃气流作用下不断被吹除，表面温度始终为相界面温度，即熔点温度，仅在固相区域存在温度分布。玻璃丝棉烧蚀过程有以下基本假设：

1)高温、高速燃气流场作用下玻璃丝棉烧蚀表面所受的燃气流气动剪切力足够大，玻璃丝棉一旦熔化即被吹除；

2)将玻璃丝棉的烧蚀过程看作是一维半无限大物体的烧蚀；

3)玻璃丝棉材料属性不受燃气流烧蚀作用的影响，始终保持不变。

采用坐标转换方法推导玻璃丝棉烧蚀速率公式，以内壁面某点为坐标原点，沿厚度方向建立的一维坐标系为定坐标系，以烧蚀界面上该点为原点建立动坐标系，初始时刻两坐标系完全重合。设烧蚀界面以速率u沿x轴正方向移动，定坐标系中的温度为T(x,t)，x轴正方向某点坐标为A(x)，随着烧蚀界面的移动，该点在动坐标系中的坐标变为A(ξ)，ξ=x-ut，则

(8)

(9)

动坐标系下玻璃丝棉烧蚀导热微分方程变换为

(10)

(11)

边界条件为

Ts(ξ,t)|ξ=0=Tp，

(12)

Ts(ξ,t)|ξ=∞=Ti.

(13)

动坐标系下玻璃丝棉烧蚀微分方程特解为

(14)

因此，玻璃丝棉固相温度场分布函数为

(15)

烧蚀界面处玻璃丝棉吸收的热量，一部分用来加热其表面，造成熔化烧蚀，一部分通过热传导方式传递到未烧蚀部分。由导热基本定律可知，热传导在玻璃丝棉烧蚀界面上热流密度为

ρscpu(Tp-Ti).

(16)

式中：λ为玻璃丝棉导热系数；cp为玻璃丝棉定压比热容。

根据能量守恒定律，烧蚀界面总的热流密度为

(17)

由(16)式和(17)式可得玻璃丝棉烧蚀速率为

(18)

因此，玻璃丝棉烧蚀深度即相界面位置可通过(19)式计算得到：

(19)

式中：tpl为玻璃丝棉内壁面温度达到材料熔点的时刻；tph为内壁面温度降到材料熔点以下的时刻。

2 玻璃丝棉消音层冲蚀损伤仿真平台

依据火箭地下井热发射燃气流场数据，利用Fortran语言编写插值程序，并通过Abaqus软件子程序开发，搭建基于Fluent、Abaqus软件的玻璃丝棉消音层冲蚀损伤仿真平台，为后续数值仿真模拟奠定基础。

2.1 燃气流场数据提取与预处理

利用Fluent软件后处理模块提取文献[2]中火箭出井过程中消音层内壁面附近燃气流场数据，包括温度、压力、单元面、节点数据等。由于文献[2]采用了域动分层动网格技术，每个时间步都会对计算域网格进行重新划分与编号，直接插值计算可能会导致消音层内壁面附近流场网格节点与单元面编号等信息发生改变。因此，为避免后续插值计算数据紊乱，采用Python语言脚本程序对数据进行提取和预处理，流程如图1所示。

图1 数据提取与预处理流程图

2.2 耦合交界面数据传递

消音层燃气冲蚀问题属于多场耦合问题，在进行数值求解时，耦合交界面的数据传递是冲蚀损伤分析的关键一步。在燃气流场数值分析时，为提高精度，网格划分比较精细，对近壁面网格进行了加密处理，如果采用相同的网格划分方案，会大大增加计算量。鉴于同精度条件下结构计算对网格数量要求远低于流场计算，所以采用不同的网格划分方案，但这样会造成耦合交界面上流体域与固体域网格密度不同、节点不匹配等问题。目前，耦合交界面不同网格间物理变量信息传递主要采用插值计算方式，这里采用反距离加权法实现插值计算。

反距离加权法是指与待插值点距离越近的点，对其产生的贡献量越大，即所有已知点对插值点产生的贡献量与其距离呈反比，可由(20)式表示：

(20)

式中：Z为插值点的估计值；Zi为第i个样本点数值；Di为样本点与待插值点的距离；p为距离的幂，相关研究结果表明，p值大小显著影响插值计算精度，p值越大，插值精度越高。本文插值计算时采用的样本数为4，幂取2.

此外，文献[2]采用Fluent软件，而本文结构场计算采用Abaqus有限元软件，这将导致二者时间步也同样存在不匹配问题。因此，采用简单的线性插值算法进行计算，估计值由(21)式得出：

(21)

式中：tl和th分别为插值时刻t的前、后燃气流场计算时刻；f(tl)和f(th)分别为tl和th时刻的样本值。

多场耦合问题数值求解时，耦合交界面数据传递方法的鲁棒性、精度关系整个数值求解的计算效率和精度，为验证所采用插值算法的可靠性，沿地下井轴向不同高度选取M、N两点，两个位置流场数值仿真结果和插值计算结果对比如图2所示。由图2可知，M、N两个位置处插值前后温度曲线均吻合较好，仅在某些时刻存在较小差异，且相对差异不大于0.56%.

图2 插值前后温度曲线

2.3 Abaqus子程序开发

虽然Abaqus分析软件提供了很多求解模型与单元库，但由于消音层烧蚀是非均匀分布载荷，所以需要单独开发用户子程序，并在任务提交时嵌入分析模型，计算时Abaqus主程序自动调用。需要开发的子程序主要包括：

1)利用FILM子程序指定出井过程中耦合交界面的对流换热系数和来流温度；

2)利用UMESHMOTION子程序定义ALE网格节点的移动速率或位移量，这里用来控制火箭出井过程中玻璃丝棉的烧蚀规律，速率由(18)式计算得出；

3)利用USDFLD子程序实现玻璃丝棉烧蚀过程中烧蚀深度、烧蚀速率等数据的可视化显示；

4)利用UEXTERNALDB子程序将井筒壁面处燃气流场所有时间步的单元信息、节点信息、温度数据等在分析开始前读入计算机内存，实现后续插值计算时流场数据的快速访问。

3 玻璃丝棉烧蚀损伤与影响因素分析

3.1 玻璃丝棉有限元模型

玻璃丝棉消音层为双层结构，外侧为普通超细玻璃棉毡，靠近燃气流一侧为防水超细玻璃纤维絮状物，外罩金属网固定，焊接在骨架上。为提高计算效率，采用子模型计算，并将其简化为平板结构，采用全六面体单元对几何模型进行结构化网格划分，靠近燃气流一侧玻璃丝棉采用等厚度网格，外侧采用渐变厚度网格，得到模型如图3所示，节点总数为65 025，单元总数为60 000，单元类型采用温度-位移耦合单元C3D8RT.

图3 玻璃丝棉网格模型及约束条件

实际安装条件下，玻璃丝棉通过焊接在骨架上扁钢固定，边界条件设置时，将其简化为背面完全固定、侧面法向支撑，如图3(b)所示。来流温度与对流换热系数大小通过调用FILM用户子程序实时获取，玻璃丝棉初始温度设定为27 ℃，玻璃丝棉材料参数如表1所示。

表1 玻璃丝棉材料参数

3.2 玻璃丝棉烧蚀损伤分析

为分析火箭出井过程中消音层燃气冲蚀损伤周向分布差异，在消音层最低端高度位置，沿周向选取子模型A1～A4. 沿轴向5个不同高度位置选取子模型B1～B5，A1与B1位置重合，示意图如图4所示，具体位置如表2所示。

表2 子模型位置分布

图4 子模型位置示意图

提取5 s时刻A1～A4位置玻璃丝棉烧蚀深度云图，如图5所示，周向不同位置玻璃丝棉烧蚀深度差异较小，烧蚀深度均在[2.38 mm,2.79 mm]。从仿真结果数据库中提取4个位置玻璃丝棉平均烧蚀深度数据，分别为2.53 mm、2.58 mm、2.52 mm、2.64 mm，相对差异不大于4.76%.

图5 周向不同位置烧蚀深度云图

周向位置不变情况下，玻璃丝棉烧蚀损伤沿地下井轴向分布特性，如图7所示为5 s时刻B1～B5位置玻璃丝棉烧蚀深度云图。从图7可以看出，沿发射井轴向玻璃丝棉烧蚀损伤差异很大，烧蚀深度分别在[2.38 mm，2.79 mm]、[2.17 mm，2.58 mm]、[1.96 mm，2.38 mm]、[0.71 mm，1.13 mm]、[0.50 mm，0.92 mm]。

图7 轴向不同位置烧蚀深度云图

提取轴向不同高度位置玻璃丝棉平均烧蚀深度数据，并绘制平均烧蚀深度随轴向高度变化曲线，如图8所示。火箭地下井发射时，消音层主要材料玻璃丝棉普遍发生熔化烧蚀，且越靠近井筒底部，玻璃丝棉烧蚀越为严重，玻璃丝棉平均烧蚀深度超过2 mm.随着轴向高度的增加，玻璃丝棉烧蚀损伤程度逐渐减小。

图8 平均烧蚀深度随轴向高度变化曲线

图9给出了轴向不同高度位置玻璃丝棉平均烧蚀速率随时间变化曲线。从图9可以看出，5条平均烧蚀速率曲线差异较大，平均烧蚀速率峰值分别为1.13 mm/s、1.07 mm/s、0.98 mm/s、1.06 mm/s、1.10 mm/s，且在燃气流场回流阶段(回流和失稳的流场特征见参考文献[2])，B1～B3位置玻璃丝棉平均烧蚀速率峰值出现的时刻逐渐推迟。此外，B1～B3位置玻璃丝棉烧蚀损伤主要发生在燃气流场回流阶段和失稳阶段，B4、B5位置玻璃丝棉烧蚀损伤则主要发生在燃气流场失稳阶段。

图9 轴向不同位置平均烧蚀速率曲线

3.3 玻璃丝棉烧蚀损伤影响因素分析

地下井建设使用过程中面临着很多实际工程问题，如不同地下井的井筒直径存在一定差异；有的增加了减振装置，导致燃气射流冲击高度增加；有的火箭发动机参数发生了改变。这些问题可能会导致井内燃气流场环境发生改变，选取烧蚀损伤最严重区域的A1位置玻璃丝棉为研究对象，开展井筒直径、冲击高度和建压速率3个因素对玻璃丝棉烧蚀损伤的影响分析，不同工况条件下各影响因素参数设置如表3所示，表中D0和H0分别为地下井原始井筒直径和冲击高度。

表3 影响因素参数设置

3.3.1 井筒直径对玻璃丝棉烧蚀损伤的影响

图10给出了不同井筒直径工况条件下火箭出井过程中的玻璃丝棉平均烧蚀速率变化曲线。

图10 不同井筒直径下平均烧蚀速率曲线

从图10中可以看出：井筒直径为D0与(D0-0.2)工况条件下，玻璃丝棉平均烧蚀速率随时间变化趋势基本一致；然而，井筒直径为(D0-1.0)工况条件下，1 s时刻前，平均烧蚀速率曲线峰值出现的时刻明显早于其他两种工况，而且平均烧蚀速率峰值也明显增大，1 s时刻以后，3个工况条件下，玻璃丝棉平均烧蚀速率曲线的幅值差异较小。此结果表明，减小井筒直径加快了燃气流场回流阶段玻璃丝棉的烧蚀速率，但对引流阶段和失稳阶段的玻璃丝棉烧蚀速率影响较小。

3.3.2 冲击高度对玻璃丝棉烧蚀损伤的影响

图11给出了不同冲击高度工况条件下，火箭出井过程中，玻璃丝棉平均烧蚀速率变化曲线。从图11中可以看出，3种工况下，玻璃丝棉平均烧蚀速率曲线随时间变化规律相似，1 s时刻前，随着冲击高度的增大，平均烧蚀速率峰值出现的时刻不断延迟，而且平均烧蚀速率峰值也明显降低。产生该现象主要是因为随着冲击高度的增大，发动机尾部喷射出的高温、高速燃气流的自由段变长，且能够得到更充分的发展，当燃气流撞击到导流锥后，被反射回来的燃气流的速率、质量率均较小，使得井内热环境逐渐趋好。但是，1 s时刻以后，3种工况下，玻璃丝棉平均烧蚀速率曲线的一致性较好。

图11 不同冲击高度下平均烧蚀速率曲线

3.3.3 建压速率对玻璃丝棉烧蚀损伤的影响

图12给出了不同建压速率工况条件下，火箭出井过程中，玻璃丝棉平均烧蚀速率变化曲线。由图12可知：建压速率为42.5 MPa/s与45.3 MPa/s工况条件下，玻璃丝棉平均烧蚀速率差异不大；但是，当建压速率增加至51.3 MPa/s时，1 s时刻前，平均烧蚀速率曲线峰值出现的时刻明显早于其他两种工况，而且平均烧蚀速率峰值也明显增大，1 s时刻以后，3个工况条件下，玻璃丝棉平均烧蚀速率曲线的幅值差异较小。此结果表明，燃气流场回流阶段，玻璃丝棉烧蚀损伤受建压速率的影响较大。

图12 不同建压速率下平均烧蚀速率曲线

4 玻璃丝棉冲蚀试验验证

为了验证玻璃丝棉冲蚀数值仿真模型的正确性，进行了玻璃丝棉的烧蚀原理性试验。试验用到的是氧-乙炔烧蚀试验台，如图13所示。焰流中心最高温度可达3 600 K左右，射流速度为50～300 m/s.

图13 氧-乙炔烧蚀试验台

为尽量模拟地下井消音层的结构，采用直径1 mm、网格间距10 mm的20号钢铁丝网对3块容重80 kg/m3、1块容重20 kg/m3玻璃丝棉块组成的20 cm厚烧蚀试件进行捆绑固定，压紧压实后试件厚度约为7 cm.将金属网固定后的试件用细铁丝绑在试验台夹具上进行烧蚀试验。仿真和真实工况中燃气流场的温度是连续变化的，而试验过程中很难使温度连续变化，实际试验时温度提前设定好，为了提高验证效果，根据流场仿真结果，将试验温度定为1 000 K、1 500 K、2 000 K，烧蚀时间定为0.5 s、1.0 s、1.5 s，共3组，每组3个试件，烧蚀比较严重的试件如图14所示。由图14可知，玻璃丝棉和金属丝网均发生熔化，熔化后的玻璃和铁丝形成了众多绿色、黄色大小不等的玻璃珠，这与历次发射试验现象一致。

图14 玻璃丝棉试件烧蚀结果

由于试验台火焰直径小，不容易模拟整体对试件进行烧蚀，所以很难用烧蚀速率表征，这里采用烧蚀前后玻璃丝棉的净质量、烧蚀量和烧蚀深度进行表述，如表4所示。

表4 玻璃丝棉烧蚀结果统计

由表4可以看出，玻璃丝棉普遍发生烧蚀，温度越高烧蚀越严重，随着烧蚀时间的增长，烧蚀量也相应增长，一定程度上验证了数值仿真模型的正确性，试件的最大烧蚀深度在2.81 mm左右，数值仿真平均烧蚀深度在2 mm左右。

5 结论

本文主要研究了地下井火箭发射过程中玻璃丝棉消音层的燃气冲蚀损伤问题，建立玻璃丝棉烧蚀损伤数学模型，搭建仿真平台并进行了数值仿真，分析不同因素对玻璃丝棉烧蚀的影响，最后进行了玻璃丝绵烧蚀原理性试验。主要得出以下结论：

1)对出井过程中玻璃丝棉烧蚀损伤特性进行了研究，结果表明：①越靠近井筒底部，玻璃丝棉受高温燃气影响越显著，底部平均烧蚀深度超过2 mm；②燃气流场温度载荷作用下，玻璃丝棉普遍发生烧蚀损伤，且主要发生在回流阶段与失稳阶段，烧蚀损伤沿井周向分布差异较小，平均烧蚀深度相对差异不大于4.76%，随着轴向高度的增加，玻璃丝棉烧蚀损伤程度逐渐减小。

2)研究了井筒直径、冲击高度与建压速率3个因素对玻璃丝棉烧蚀损伤的影响，结果表明，减小井筒直径或增大建压速率将加剧玻璃丝棉的烧蚀损伤，而增大冲击高度将减小玻璃丝棉的烧蚀损伤。

3)进行了玻璃丝棉的冲蚀原理试验，结果表明在燃气流的作用下，玻璃丝棉普遍发生烧蚀，说明了仿真模型的正确性，对后期消音层的处理及修复有一定指导意义。