水下固体火箭发动机推力脉动特征研究

唐云龙，李世鹏，刘　筑，隋　欣，王宁飞

(北京理工大学 宇航学院，北京　100081)



水下固体火箭发动机推力脉动特征研究

唐云龙，李世鹏，刘筑，隋欣，王宁飞

(北京理工大学 宇航学院，北京100081)

对水下固体火箭发动机的动态特性进行了研究，根据动网格及VOF(Volume of Fluid)基本原理，构建了用于水下固体火箭发动机动态射流响应的轴对称模型；根据壁面积分基本原理及被动方式求解水下工作发动机推力的计算方法，成功模拟了发动机的受力运动；对不同工况下发动机工作情况进行仿真，研究了水下发动机动态模型下的射流结构，发现在来流作用下，射流气体会受到来流边界层的挤压，产生颈缩现象，之后逐渐向外膨胀；研究结果显示，在来流边界层作用下，射流所产生的激波在向下游发展之前，将经历多次折射，该研究结果与文献中相关实验结果较吻合；对发动机不同工况下的推力脉冲峰值进行分析，找出了动态模型下发动机点火初期推力脉冲峰值特征随工作环境的变化规律。最后，对有/无轴向浮力时的发动机动态特征进行了对比分析。

动网格；水下射流；边界层；推力峰值

0　引言

利用动网格方法，研究水下固体火箭发动机运动过程的射流结构及推力问题，与水下工作发动机的真实情况更加接近，且有利于掌握水下火箭的工作规律，能为水下发动机设计及工程应用提供帮助。国内外在水下射流流场结构与流动特性研究领域发展均较缓慢[1]。唐嘉宁等[2]将球形气泡模型与CFD计算结果进行比较，发现球形气泡模型结构简单、计算速度快,并能有效反映发动机点火瞬间的推力峰值特性。但该模型基于气体在整个射流过程中以单一气泡结构存在，并沿径向等速无限膨胀假设，并未考虑水-气间相互作用，因此对推力的解释并不是很全面。Lindau等[3]对考虑空泡作用下高速超空化水下推进过程进行了数值研究，提出了燃气、水和水蒸气的三相流计算模型，但未对水下推进过程中的流动特性及发动机工作性能做进一步分析。唐嘉宁等[4]建立了水下超音速气体射流轴对称数学模型，并将数值计算结果与实验结果进行对比，验证了该模型的合理性和准确性[5]。

为研究固体火箭发动机在水下工作时水的压差阻力和来流速度等对其稳态推力的影响，本文建立了发动机水下流场模型,并在发动机前端加设弹体段，通过数值方法，分析了固体火箭发动机水下运动达到平衡速度后的流场结构及发射初期的推力振荡峰问题。研究结果对导弹水下发射的流体动力特性分析及发动机设计具有指导意义。

1　基本模型及方程

1.1VOF原理

VOF是通过研究网格单元中流体和网格体积比函数F对自由面进行确定，并追踪流体的变化，并非追踪自由液面上质点的运动。VOF方法可处理自由面重入等强非线性现象,其有计算时间短、存储量少等特点，但在处理F变化时稍显繁琐，其有一定人为因素。若F=1，则该单元为指定相流体单元；若F=0，则该单元为无指定相流体单元；当0

1.2动网格技术

动网格技术方便研究弹体与水间相对运动，又能保证网格块之间流动变量的插值精度。在使用动网格技术时，一般将流动域划分为2个或2个以上的计算域，相邻的计算域由1个或多个交界面连接，并沿交界面滑动。计算过程中，运动域沿着交界面滑动，每一时间步后，重新计算网格位置，确定交界面上网格点的相对位置，并在交界面处进行通量传递[7]。

1.3层变网格方法

在四边形网格区域内，动态层变方法可根据与运动相邻的网格层高度增加或减少网格层数。在计算中，动网格模型可指定一理想的高度。如图1所示，邻近边界的网格单元层根据层j的单元层高度h来分裂或与邻层合并成新层。若层j中单元体积处于膨胀状态，其膨胀极限如下[8]：

(1)

式中hideal为理想单元高度；ah为全局单元层分裂或合并因子。

当h>(1+ah)hideal时，层j中的单元将分裂成高度hideal的单元层和高度为h-hideal的单元层。若层j中单元体积处于被压缩状态，则当压缩到h

图1　动态网格层

1.4动网格控制方程

φ代表通用标量，静止区域和运动区域内的方程可写为如下形式积分方程[9]：

(2)

式中ρ为密度；ug为网格运动速度；Γ为湍流扩散率；Sφ为源项；φ为运动区域内的任意控制体；∂V代表控制体的边界。

方程(2)的时间导数项可利用一阶向后差分写成如下形式：

(3)

其中，n及n+1表示时间。第n+1时间的层网格体积可通过式(4)计算：

(4)

其中，dV/dt为控制体的体积时间导数。为满足网格守恒率，控制体的体积时间导数通过式(5)进行计算：

(5)

(6)

其中，∂Vj为整个时间步Δt上因控制容积面j膨胀而引起的体积变化。

2　计算方法

2.1推力

在计算过程中，由于动网格边界运动规律是未知的，故利用被动方式来求解水下工作发动机的推力。其中，Fi和Fe分别表示作用在发动机外表面和内表面上压强的合力(如图2所示)，则推力表达式如下[10]：

(6)

图2　发动机计算方法

2.2速度

本文研究的水下发动机涉及被迫运动问题，其受力非常复杂。研究过程中，利用轴对称模型，即发动机只仅在轴线方向运动。根据牛顿定理可知轴线上力F可表示为：

F=ma

(7)

故某一时刻速度表达式为：

(8)

式中v(NΔt)为NΔt时刻的速度；Δt为计算时间步长；Fk为kΔt时刻所对应的推力；M为发动机总质量；Δm为推进剂在单个时间步长内的质量损失。

3　仿真模型

在结合水下发动机的实际工作情况的基础上[9]，建立了相应的二维轴对称模型。计算域如图3(a)所示，包含喷管内流区域和外部射流区域。发动机喷管内流区域及邻近发动机外流区域使用非结构网格作为过渡。在远离发动机区域使用结构网格作为滑动区域。

(a)计算网格

(b)发动机周边网格　　　　(c)前端不动网格

图3(a)中，上下的滑动网格与静态网格采用接触边界进行接触。在运动区域的起始与结束段采用内部网格作为动态和静态网格的分界(如图3(c)所示)。整个区域的外边界均为压力出口，而发动机喷管入口则采用压力入口作为边界条件。因为采用了滑动网格技术，网格的产生、消失及变形都会对发动机进流场产生较大的影响。所以，为了减少动网格本身对流场的影响，动网格与静网格的接触区域选在流场远端，且满足长时间非稳态计算的要求。

4　仿真结果及分析

4.1流场

4.1.1流场结构特征分析

图4为利用图3中的水下发动机单自由度运动模型达到平衡速度时的计算结果及发动机水洞实验结果。两者都是采用燃烧室内压恒定的射流模型，且在来流作用下，轴线推力达到平衡值，都趋近于零。

(a)水下发动机射流体积分数仿真结果

(b)水下发动机射流体积分数实验结果[11]

(c)水下发动机射流密度分布仿真结果

(d)水下发动机射流激波分布仿真结果

从来流与水下发动机射流耦合作用下的射流及激波结构上分析，可知图4 (a)、(b) 中射流仿真的体积分数分布与文献[11]中仿真结果(如图4(c)所示)及文献[12]中实验结果较一致。由于受到来流惯性作用，喷管出口处会产生边界层，射流形成的激波(如图4(d))在边界层内不断折射导致通道变窄(如图中A、B、E、F所示)气体膨胀受阻随着射流的方向发展。来流对射流的作用减弱，射流开始向外膨胀(如图中C、D、G、H所示)，其中A、B、E、F为燃气射流起始段(如图(a)中L1所示)，C、D、G、H之后为燃气射流主体段(如图(a)中L3所示)，在两者之间的为过渡段(如图(a)中L2所示)。通过模拟水下射流流场的运动过程，发现在有来流情况下，由于外界水流场作用，燃气射流不再断裂，并未发生回击现象，与文献[11]中结论一致，证明该模型能真实反映动态过程。

图5为射流流动特性与轴线上压强和速度无因次化后的对应关系。按图4(a)中的方法对轴线上的压强和速度进行划分，研究发现起始段、过段段和主体段特征较明显，颈缩对射流进行了二次加速。

4.1.2水下发动机运动平衡后轴线特征

何淼生[13]、王春旭[14]、曹嘉怡[15]等学者已经对水下发动机的开机、关机及出水过程进行了详尽研究。所以，本文将不再对其进行重复性研究，只对在水下发动机运动平衡后的特征进行非稳态分析。

图6为水下发动机在达到平衡速度后，在不同时刻发动机近流场的分布特征。如图6(a)、(b)所示，水下发动机射流流场在边界层与激波共同作用下而产生的颈缩现象。

图5　某一时刻射流流动特性与轴线上压强和速度无因次化后的对应关系

(a)压强分布

(b)马赫数分布图

从图6(a)、(b)可知，当发动机速度达到平衡后，尾流场结构及发动机前段流场均无明显改变，射流起始段长度基本固定。采用图5中方法结合流场变化，可描述不同工况下发动机尾流场中颈缩点与射流起始段长度，因为衡后射流结构及发动机轴线特征不再变化，所以不再对发动机速度平衡后时刻的流场分布进行特殊说明。

4.2推力特征分析

图7(a)为固体火箭发动机推力-时间曲线仿真结果，如图7所示，其动态模型的推力无大幅振荡，与静态模型有明显差别。在射流初期，射流喷出喷管形成气泡，并产生强劲的推力脉冲，本文仿真模型与文献[13]中实验模型均采用端燃药柱，认为燃烧室内压不变，且对径向位移添加了相应约束，只在轴向上单一自由度运动。

(a)航行器推力-时间曲线

(b)航行器粘性阻力-时间曲线

(c)航行器速度-时间曲线

(d)航行器位移时间曲线

通过对水下发动机轴向推力的监测数据比较，发现本文的仿真结果与文献[13]中实验结果一致；随着速度和推力的增加，发动机在压差阻力作用下，开始减小直至到零，并在零值附近做小幅值振荡。如图7(b)所示，发动机的黏性阻力相对其压差阻力可忽略不计。如图7(c)所示，发动机运动速度在推力脉冲的作用下不断增加，当来流阻力与发动机推力平衡时，速度也趋于平衡。如图7(d)所示，通过速度对时间积分，可得发动机位移变化曲线，当发动机达到平衡速度后，位移曲线基本呈线性改变。

为了研究水下固体发动机单自由度运动模型在不同工况下发射初期推力的振荡问题，本文将不同总压、总温、水深、扩张比及质量按照单调变化趋势分成5种工况进行比较(见图8)。从图8(a)可知，水下固体火箭发动机在点火初期会出现一个非常大的峰值，该峰值对水下发动机的运动特征具有很大影响。

从图8(a)～(e)可知，在计算工况下，总压越大，推力脉冲峰值越大；水深越大，推力脉冲峰值越小；扩张比越小，推力脉冲峰值越大；总温越高，推力脉冲峰值越大；质量越大，发动机推力脉冲峰值越大，且推力脉冲与发动机推力积分及平衡速度之间具有相关性。

不同质量下的速度-时间曲线见图9。由图9可知，水下发动机质量变化不会改变其运动的平衡速度，但在运动相同时间下，水下固体发动机的位移却随质量的增加而不断减小，这说明水下发动机运动速度的爬升过程不同。速度爬升段的面积增加与质量减少(如图10所示)基本呈线性关系，且是否考虑浮力对其增量变化率基本没有影响。

(a)燃烧室压强对发动推力峰值的影响　　　　　　　(b)水深对发动机推力峰值的影响

(e)扩张比对推力峰值的影响

图9　不同质量下的速度-时间曲线　

图10　质量减少量与位移增量的关系

5　结论

对发动机在不同工况下的运动过程进行了仿真研究，发现水下发动机在运动过程中来流会与发动机尾部燃气射流发生耦合；在发动机运动模型中，尾流不会产生断裂/回击；来流速度与射流激波的相互作用会在发动机尾部形成边界层，使射流形成颈缩，且射流进行二次加速；根据射流特征对燃气射流的起始段、过渡段和主段进行划分，并与轴线上的速度和压强特征进行对应；通过对不同时刻轴线及速度特性分析，发现当发动机达到平衡速度后，其运动规律及射流特性基本保持不变；动态模型发动机推力在初期会形成一个脉冲，且粘性阻力相对推力可忽略不计；通过对不同工况下发动机推力特征分析，得出以下结论：水深越大，发动机推力脉冲越小；燃烧室压强越大，发动机推力脉冲越大；温度越大，发动机推力脉冲越小；扩张比越大，发动机推力脉冲越大。改变发动机的质量不会改变发动机的平衡速度，但发动机质量越小，加速越快，且发动机轴向浮力对其速度增量影响很小。

[1]甘晓松，贾有军.水下燃气射流流场数值研究[J].固体火箭技术，2009，31(1):23-26.

[2]Tang J N,Li S P,Wang N F,et a1.Flow structures of gaseous jet injected into liquid for underwater propulsion[C]//46th AIAA/ASME/ SAE/ASEE Joint Propulsion Conference and Exhibit,Nashville,IN.2010,AIAA 2010-6911.

[3]Lindau J W,Venkateswaran S,Kunz R F,et a1.Multiphase computations for underwater propulsive flows[C]//16th AIAA Computational Fluid Dynamic Conference,Orland,Florida,AIAA 2003-4105.

[4]Tang Jia-ning,Wang Ning-fei,Wei Shyy.Flow structures of gaseous jets injected into water for underwater propulsion[J].Acta Mechanica Sinica,2011,27(4):461-472.

[5]唐嘉宁，李世鹏，王宁飞.水下固体火箭发动机的负推力现象研究[J].固体火箭技术，2012，35(3)：325-344.

[6]康玉东，孙冰．燃气非平衡流再生冷却流动传热数值模拟[J].推进技术，2011，32(1):119-124.

[7]粱伟,宋文萍.基于滑动网格的螺旋桨复杂绕流数值模拟研究[J].航空计算技术,2011,42(3):45-49.

[8]杨猛.计算流体力学两相流流动的模拟及两相流模型的研究[D].天津大学,2005.

[9]杨晓光,陈焕龙,刘华坪,等.导弹水下运动及出水过程的三维流场仿真[J].弹道学报,2010,22(1):107-110.

[10]唐嘉宁,刘向阳.水下固体火箭发动机推力特性研究[J].导弹与航行运载技术,2012(5):15-21.

[11]何晓.喷气推进导弹空泡流的尺度效应研究[D].上海交通大学,2009.

[12]Ivan N Kirschner,Neal E Fine,James S Uhlman,et a1.Supercavitation research and development[J].Undersea Defense Technologies,Waikiki,HI,2001:1-10.

[13]何淼生,覃粒子,何佳磊,等.水下超声速喷嘴起动/关机过程射流的动态不稳定性研究[J].推进技术,2014 ,35(4):523-529.

[14]王春旭,张涛,侯国祥.基于涡声理论的水下射流噪声预报[J].船舶力学,2010，14(6):670-677.

[15]曹嘉怡,鲁传敬,李杰,等.潜射导弹水下垂直自抛发射过程研究[J].水动力学研究与进展,2006,21(6):752-759.

(编辑：崔贤彬)

Research on thrust fluctuation characteristics of underwater solid rocket motor

TANG Yun-long，LI Shi-peng，LIU Zhu，SUI Xin，WANG Ning-fei

(School of Aerospace Engineering,Beijing Institute of Technology,Beijing100081,China)

The dynamic property of underwater solid rocket motor was studied in this article.Utilizing dynamic grid and the VOF principle,an axisymmetric model describing the dynamic response of the fluidic thrust of underwater solid rocket was established.Movement of the motor was successfully simulated based on the principle of surface integral and the computational method of underwater motor's thrust.By simulating motors under different cases using the dynamic model,the jet structure of underwater motor was studied.It turn out that the jet gas is extruded by the flow in the boundary layer,and generated necking phenomenon,and then expanded outward gradually.The results indicate that the shock wave produced by the jet will refract several times before going downstream,which is consistent with experimental results of other literature.The thrust's pulse peaks in different cases,were analyzed,and the change law of the features of underwater motors with the variation of environment was obtained.Finally,the dynamic characteristics of motors in conditions with/without buoyancy were compared and analyzed.

dynamic grid;water jet;boundary layer;thrust peak

2015-06-04；

2016-01-18。

国家自然科学基金(11072032)。

唐云龙(1986—)，男，博士生，从事固体火箭发动机仿真及测试研究。E-mail：ks669@126.com

V438

A

1006-2793(2016)04-0476-06

10.7673/j.issn.1006-2793.2016.04.005