冲击作用下推进剂变形的流固耦合分析方法

金铭君，李 强

(西北工业大学 燃烧、热结构与内流场重点实验室，西安 710072)

冲击作用下推进剂变形的流固耦合分析方法

金铭君，李 强

(西北工业大学 燃烧、热结构与内流场重点实验室，西安 710072)

固体火箭发动机的点火过程是一个复杂多变的理化过程，具有时间短、升温、升压梯度大等特点。针对固体火箭发动机点火过程中的装药结构完整性问题，文中建立了一套用于分析冲击作用下固体推进剂变形现象的仿真模型。采用RANS和ALE方法，分别对流体域和固体域进行求解，以两场独立交叉耦合迭代的模式实现了仿真过程。以一个推进剂冷流冲击实验作为算例，对仿真模型进行了验证，计算值与测量值间误差不超过10%，仿真模型计算可靠，具有向固体火箭发动机实际点火过程拓展的价值。

固体火箭发动机；点火过程；推进剂变形；流固耦合

0 引言

固体火箭发动机的点火过程是一个短暂的快速升压过程。其中，高温高压燃气对固体推进剂装药的冲击作用，会导致装药结构的变形、破损和移位，甚至导致发动机失效爆炸。因此，为确保整个发动机的工作性能及可靠性、安全性，对发动机点火瞬态过程开展装药结构完整性及瞬态流畅变化研究是非常有必要的。一直以来，固体火箭发动机的研制和定型主要以试验为主，周期长、耗资大，加之燃烧室内工作环境恶劣，点火过程又极为短暂，且复杂多变，对其进行实验研究是极为困难的。近年来，随着各种高精度、高效率计算方法的提出和完善，以及高速计算机的快速发展，计算机数值仿真技术凭借其高效、经济、安全等优势，在火箭发动机的研制中，占据了愈发重要的地位。

国内外大量研究表明，升压梯度是破坏装药结构完整性的重要因素[1]。在点火阶段，过高的压强峰或升压速率会对推进剂药柱产生冲击，造成药柱变形或移位，而药柱变形或移位又使得流动区域发生了变化。对于这样一个典型的流固交互问题，采用流固耦合法对其分析，无疑是最佳选择。1982年，Tante G L将流体和结构之间的耦合关系作为动力学模型，验证了燃烧室和药柱的相对运动，建立了用于分析SRM点火爆炸原因的二维非定常模型。1991年，Titan IV SRMU地面试车失败，掀起了针对点火冲击影响的流固耦合方法研究热潮。1991年，Johnston采用二维流场耦合结构变形模型，计算了内压作用下药柱变形及其对流场的影响，给出了较为安全的药柱临界模数[2]。此后，许多研究小组开展了此方面的研究工作，并得到了有价值的结论。国内报道的主要有海军航空工程学院飞行器工程系通过MPCCI连接FLUENT和ABAQUS软件，对SRM快速升压过程的研究[1]；以及南京理工大学机械工程学院利用CFX和ANSYS软件，对固体推进剂裂纹内点火过程的流固耦合研究等[3]。

本文着眼于具有良好适应性的高精度SRM点火过程流固耦合方法研究，提出了一种适用的流固耦合计算分析模型。由于实际发动机点火过程的推进剂应变和变形数据无法通过实验获得，因此本文以一个推进剂冷流冲击实验作为算例，对计算模型进行了校验，通过对比计算结果与实验结果中测试段内的压强曲线和推进剂变形最大位移的大小，来对仿真模型的准确性做出评价。

1 计算模型

1.1 流动计算模型

针对气流对推进剂冲击过程的流动特性，对流动模型做如下假设：

(1)计算中所涉及流体为有粘可压缩牛顿流体，且可视为各向同性的连续介质；

(2)冷流冲击过程中流体温度脉动较小，忽略由温度脉动带来的粘性系数的脉动；

(3)流场运算中忽略重力影响。

将密度ρ、速度ui、压力p及温度T经上式处理后，代入有粘可压缩N-S方程组，即可得到可压缩雷诺控制方程组，进行必要简化后，有[8]

(1)

(2)

(3)

1.2 湍流模型

本文采用涡粘性方法中的一方程模型：Spalart-Allmaras模型(S-A模型)来求解雷诺平均方程。

对于一方程模型的一般形式，有

(4)

式中SP为生成项；SD为耗散项；D为扩散项。

(5)

其中

这里，d为到物面的最近距离。对于式中其他参变量的关系，有

g=r+Cw2(r6-r)

上述各式中，常系数取值如表1所示。

表1 S-A模型中常系数取值

1.3 推进剂结构力学计算模型

对推进剂结构模型做如下假设：

(1)认为冲击时间短暂，固体推进剂的粘弹特性不能完全体现，将其近似为大变形弹性模型；

(2)将推进剂材料视为均质的各向同性Kirchhoff材料，忽略重力影响。

对于上述假设，材料本构关系[9-10]满足：

Sij=Cijk1Ek1=λEkkδij+2μEij

(6)

其中，Sij为基尔霍夫应力；Eij为格林应变；Cijk1为弹性模量的四阶张量；λ和μ为拉梅(Lamé)常数，可表示为

λ=υE/(1+υ)(1-2υ)

μ=E/2(1+υ)

1.4 非线性有限元方法

采用ALE(Arbitrary Lagrangian-Eulerian)有限元方法，对结构域进行离散。在ALE方法中，由于质量守恒被强加作为偏微分方程，所以需要建立一种弱形式[10]。将ALE域划分为单元，定义所有非独立变量为单元坐标的函数，对于单元e，其ALE坐标为

χ(ξe)=Φe(ξe)=χ1N1(ξe)

其中，ξe为单元e的坐标。对于网格运动，有

其中，x1(t)为节点的运动。对于网格速度，有

对于有限元矩阵方程，有连续方程：

其中,Mρ、Lρ、Kρ分别为容量、转换和散度矩阵：

对于动量方程，有：

其中M和L分别是广义质量和传递矩阵，对应于在参考构型描述下的速度fint和fext分别为内力和外力向量：

1.5 耦合计算

目前，流固耦合问题的求解方式有两类：两场交叉迭代和直接全场同时求解。后者计算量巨大，且具体实现十分复杂，故本文采用两场分开进行交叉迭代计算，各自单独使用独立解算器和网格进行解算的方式。在两场数据交换的处理上，一般是依靠在两套非匹配网格间进行双向插值来完成的，其精度受两插值网格间重合度等诸多因素影响。本文采用在交界面处通过插值重构一层虚网格进行数据传递和相关变量解算的耦合方法，数据通过最小二乘法在两场和虚拟网格间进行传递，减少了来自双向插值带来的误差；交界面节点力可在虚网格上直接进行计算，节省了计算时间和计算量。设上角标nf为流体解算器中的时间步，ns则为结构解算器中的时间步；下标f代表流体解算器中的变量，s代表结构解算器中的变量，f-s则代表交界面处的变量。

通常固体膨胀波速cd比流体中波速c大，出于对算法稳定性的考虑，固体解算器中的时间步长ts通常小于流体解算器中的时间步长tf。因此，两场求解器中分别取其最佳时间步长进行交叉迭代，其具体流程和数据传递如图1和图2所示。

2 实验

2.1 实验方法

该实验采用一套冷流来流模拟装置来模拟发动机点火过程中压强迅速升高的燃气流，对缩比的推进剂进行冲击实验。实验开始前，先在高速摄影仪驱动程序界面上对推进剂试件取一系列距离已知的标定点；实验时，通过高速摄影仪来获得实验过程的影像资料，并在测试段首尾两端配置压力传感器来对内腔压强进行监测；实验结束后，对影像资料中标定点间的距离变化情况进行分析，获得推进剂形变、形变速度随时间变化的曲线[11-13]。

2.2 实验装置

实验装置主要由3部分组成：氮气恒压供气系统，推进剂变形测量段及数据采集系统。氮气恒压供气系统由4个氮气气瓶和稳压装置组成，可提供2～4.5 MPa的冲击压强。实验测量段内腔尺寸为240 mm×200 mm×140 mm，内有固定槽，可保证在测量过程中试件没有移动；前端为30 mm厚的有机玻璃，在保证测量段密闭性的同时，又兼顾了图像采集系统对视野清晰度的要求；实验采用两路测压，测压孔配置在测量段后盖板上，测量点P1临近气流入口处，测量点P2临近气流出口位置。数据采集系统分为图像采集系统和压强测量系统两部分。图像采集部分采用美国Phanton高速摄影仪(V4.3)，所采用的传感器分辨率为800×600，采样率为300帧；压强测量部分采用Dewetron并行数据采集系统。

2.3 实验方案

实验对不同冲击压强下推进剂的变形情况进行了测量，为了减小实验误差，对每个工况进行了3次实验，最终结果为3次实验结果的平均值。本文采用前文所述计算模型，对其中冲击压强为3.3、3.6、4.1 MPa的工况开展了数值计算工作。

3 结果对比

3.1 实验结果

图4给出了在冲击压强为3.3、3.6、4.1 MPa的工况下，气流入口和出口处压强随时间变化曲线。其中，曲线p1代表气流入口处压强，p2代表气流出口处压强。

在各个工况下，冲击压强达到最大值需要约0.12 s，与实际发动机点火过程的压强建立时间大致相当，可认为该冷流推进剂冲击实验具备一定对真实发动机点火过程的参照性。

在各个工况下，气流入口与出口处存在明显的压强差，分析认为这是造成推进剂变形的主要原因之一。推进剂形变的形状与氮气冲击的位置有关，在该实验中，冲击位置位于推进剂的右侧，推进剂受冲击后向左凹陷，位移最大处为气流冲击位置的中央处，如图5所示。各工况下入-出口压强差、推进剂最大位移及形变速度如表2所示。

冲击压强/MPa前后压差/MPa推进剂形变大小/mm推进剂形变速率/(mm/s)3．30．1610．70846．9363．60．2360．81753．1304．10．3291．19259．792

3.2 计算结果

参照上述实验，依据其设计参数配置物理模型和边界条件，开展了对推进剂冷流冲击的数值计算工作。流动求解器采用了有限体积方法，空间对流项的离散采用了HLLC格式，粘性项差分采用WENO格式，时间项上采用四阶Runge-Kutta积分方法。在计算中，参照实验测量到的推进剂位移最大位置，设置了位移监测点。由于实验结果分析中只给出了氮气冲击造成的推进剂变形的最大位移结果，而未能给出具体位移情况随时间变化的过程，考虑到氮气的冲击作用只是一个脉冲过程，因此在数值计算中，在气流入口位置给定一个单脉冲过程的冲击压力作为初始边界条件。计算域分为16个分区，流体计算网格总数约为28万，其计算结果如下：

首先，由表3可知，各工况下最大位移的计算值和测量值间误差不超过10%，可认为计算结果较为准确。其次，在图6中给出的推进剂最大位移随时间变化曲线可看到，在冲击过后最初的0.02 s内，3个工况下推进剂均没有发生明显的形变。此后，推进剂开始变形，且形变量逐渐变大。对比图4可知，在冲击最初的0.02 s内，虽然氮气气源和实验器内的绝对压差很大，但实验器内气流入口和出口间的相对压差却较小，因此推进剂的变形很不明显；而在0.02 s之后，随气流入口和出口间相对压差的逐渐增大，推进剂变形逐渐明显，且在压差达到最大值的时刻，推进剂的位移变化几乎同时达到最大值。最后，在冲击过后的降压阶段，由于计算边界无法准确模拟氮气压强衰减的过程，可认为此后的计算结果没有意义。

冲击压强/MPa推进剂最大位移/mm计算值实验值3．30．7210．7083．60．8920．8174．11．0861．192

由于实验条件所限，实验结果并没有给出推进剂变形大小随冲击压强渐变的具体过程，因此在计算结果中，只有推进剂形变的最大值和实验结果具有比照作用，但其得到的过程值可认为是具有一定预验价值的。通过对比计算结果与实验结果中3个工况下的推进剂变形最大位移以及测试段内升压过程的压强曲线，验证了仿真模型的可行性和准确性。

4 结束语

(1)为了对气流冲击作用下推进剂的变形现象开展数值仿真研究，本文提出了一套基于RANS和ALE方法的流固耦合仿真模型。

(2)通过在交界面构造虚拟网格的方法，实现了两场间的数据交互，既减少了来自双向插值带来的误差，又节省了计算时间和计算量。

(3)以推进剂冷流冲击实验作为算例，对仿真模型进行了校验，计算结果与实验结果误差不超过10%，验证了仿真模型的可靠性和准确性。

[1] 于胜春, 赵汝岩, 许涛, 等.固体火箭发动机快速升压过程的流固耦合分析[J].固体火箭技术, 2008, 31(3):232-235.

[2] Johnston W A. Flow-structural analysis of the Ariane-5 solid rocket motor during ignition transient[R]. AIAA 2011-4056.

[3] 韩波, 周长省, 陈雄.固体推进剂裂纹内点火过程流固耦合数值仿真[J].固体火箭技术, 2011, 34(2):180-183.

[4] 阎超.计算流体力学方法及应用[M].北京: 北京航空航天大学出版社, 2006.

[5] Mark Salita. Modern SRM ignition transient modeling(part 1): introduction and physical models[R]. AIAA 2001-3443.

[6] 王新月, 胡春波, 张堃元, 等.气体动力学基础[M].西安: 西北工业大学出版社, 2006.

[7] David C Wilcox. Turbulence modeling for CFD[M]. California: DCW Industries, Inc., 1993.

[8] 陈矛章.粘性流体动力学基础[M].北京: 高等教育出版社, 1993.

[9] Toro E F. Riemann solvers and numerical methods for fluid dynamics: a practical introduction[M]. Berlin:Springer-Verlag Berlin and Heidelberg GmbH & Co. K, 1999.

[10] Ted Belytschko, Wing Kam Liu, Brian Moran. Nonlinear finite elements for continua and structures[M]. London: John Wiley & Sons Ltd., 1998.

[11] Soares D Jr, Rodrigues G G, Gonc alves K A. An efficient multi-time-step implicit-explicit method to analyze solid-fluid coupled systems discretized by unconditionally stable time-domain finite element procedures[J]. Computers and Structures, 2010, 88(5-6):387-398.

[12] Luciano Garelli, Rodrigo R Paz, Mario A Storti. Fluid-structure interaction study of the start-up of a rocket engine nozzle[J]. Computers & Fluids, 2010, 39(7):1208-1218.

[13] Guerraiche Bilal. Invstigation on mechanical behavior of composite propellant during ignition process[D].西安: 西北工业大学, 2013.

(编辑：崔贤彬)

A solid-fluid interaction model for propellant deformation under impact condition

JIN Ming-jun，LI Qiang

(Science and Technology on Combustion, Internal Flow and Thermal-Structure Laboratory,Northwestern Polytechnical University, Xi'an 710072, China)

The ignition of SRM is a complicated and changeful physicochemical process, which has the features of extremely short duration and particularly high temperature & pressure gradient. An applicable solid-fluid interaction simulation model was established to analyze the phenomenon of propellant deformation under impact effect for the purpose of dealing with grain structure integrity problem during SRM ignition transient. RANS and ALE method were used independently in fluid and solid sub-domain. The interaction between two sub-domains was accomplished by the transmission of interface forces. The model was validated with a cold flow test. The results of simulation agree well with test data. It shows that the model has a relatively high accuracy and could be used to investigate the ignition process in a practical SRM.

solid rocket motor；ignition transient；propellant deformation；solid-fluid interaction

2016-02-25；

2016-10-27。

金铭君(1988—)，男，硕士生，从事固体火箭发动机点火过程流固耦合研究。E-mail:jackking413@163.com

V512

A

1006-2793(2017)02-0158-06

10.7673/j.issn.1006-2793.2017.02.005