嵌银丝端燃装药耦合传热的内弹道数值分析

肖志平，韩　波，张　亮，职世君，邢国强

(中国空空导弹研究院，洛阳　471000)



嵌银丝端燃装药耦合传热的内弹道数值分析

肖志平，韩波，张亮，职世君，邢国强

(中国空空导弹研究院，洛阳471000)

针对内埋金属丝端燃装药燃烧过程中出现压强峰值的问题，耦合金属丝区域的传热过程，建立了内埋金属丝端燃装药燃速计算模型。通过对有限元软件，进行二次开发实现了内埋金属丝装药传热及燃面退移过程的模拟。结果显示，沿金属丝被加热推进剂的温度、燃速、增速比受燃烧室压强影响较大。在此基础上，进行了嵌单根银丝端燃装药内弹道仿真，得到了压强，燃面面积，沿金属丝燃速随时间变化规律。结果表明，计算模型可有效模拟内埋金属丝端燃装药燃烧过程，与实验数据吻合较好，从理论上解释了燃烧过程中出现压强峰值的实验现象。

嵌金属丝装药；二次开发；内弹道；传热；燃速

0　引言

端燃药柱通过内埋金属丝，可增加燃烧面积，达到提高总压的目的。该过程的基本原理是在装药燃烧过程中，金属丝受到高温燃气作用，温度急剧升高，从而对金属丝附近的装药进行了加热，提高了装药的温度。金属丝附近药柱温度提高，导致其燃速较其他各处要快，在金属丝处将形成锥面，随着锥面的增大，燃烧面积也不断增大，整个燃烧室的压强也相应提高。

试验显示，嵌金属丝后，由于燃面面积增大，压强会显著增加，并在初始时刻出现压强峰[1]，之后压强达到平衡。

国外King M 等[2]用数值方法较系统地研究了金属丝的材质、直径以及装药的性质、间隙等对嵌金属丝端燃药柱燃速的影响。Gossant B等[3]计算了嵌金属丝端燃装药的增速比。国内张有为等探讨了降低嵌金属丝端燃药柱固体火箭发动机点火初始推力峰值的方法[4]，并先后进行一系列嵌金属丝端燃装药的特性研究[5-6]。关于该方面的数值计算，研究人员都开展了大量仿真工作[7-9]。但由于问题的复杂性，现阶段对嵌金属丝端燃装药的燃面变化过程还缺乏深入的研究，上述文献多采用简化模型，将燃面等效为母线为直线的锥面，计算得到的内弹道曲线并没有捕捉到初始压强峰值，与试验不符。而实际过程中，燃面形貌在不断发生变化，压强会出现较高峰值，最后达到稳态时，燃烧锥面的母线才是直线。

本文通过对有限元软件进行二次开发耦合嵌金属丝端燃装药的传热和燃面变化过程，通过数值仿真，得到了压强-时间变化曲线，从理论上解释了初始压强峰值出现的原因。

1　数学模型

内埋银丝端燃装药内弹道过程中，存在着内弹道微分方程与传热方程的相互耦合，计算模型示意图见图1。

图1　嵌银丝装药模型示意图

1.1传热方程

热传导微分方程[10-11]表达形式为

(1)

式中T、ρ、c、ki、Q分别表示温度、密度、比热容、热导率、单位体积热生成率。

燃气对金属丝的热传导边界条件采用对流换热边界，热流密度[2]采用下式：

(2)

式中ρprop、rw、cp、Tgas、Tw分别为推进剂密度、沿金属丝推进剂燃速、推进剂定压比热容、火焰温度、金属丝表面温度。

在计算过程中，随着装药燃面的不断推移；包裹金属丝的推进剂燃烧完后，金属丝会暴露出来，从而形成新的与燃气相互作用的边界；另外，还要考虑金属丝在燃气中温度达到熔点后的熔化过程。

1.2内弹道计算方程

内弹道计算微分方程[12]为

(3)

式中pc、t、R、T0、Vc、S、r、ρp、CD、At分别表示燃烧室总压、时间、气体常数、总温、燃烧室自由容积、燃面面积、燃速、装药密度、流量系数、喷管喉部截面积。

装药燃速计算公式：

(4)

燃速系数a(T)采用King M提出的改进模型进行计算[2]。

(5)

式中T为推进剂温度。

金属丝对推进剂进行加热，相当于提高了附近推进剂初温T，燃速系数随推进剂温度是单调递增的，温度越高，燃速系数越大，燃速也相应的增大。通过实验可确定式(5)中的相关系数。

由于燃速随压强和装药温度的变化规律可通过实验来确定，因此不需计算火焰对装药的传热过程，避免了引入更多的不确定因素。

实验表明，增速比随金属丝直径先增大、后减小，即有获得最大增速效果的直径[2]。仿真计算表明，若直接取燃面处金属丝与装药接触点上的温度，或者接触面轴向增量燃烧厚度上的平均温度来计算燃速系数，会得到金属丝越细、增速比(沿金属丝被加热装药燃速和基础燃速比值)越大的结论，不符合试验规律。通过实验和仿真交互验算 ，本文修正了King M给出的公式[2]，采用下式作为金属丝附近装药温度：

(6)

式中α为推进剂的热扩散系数；rwire为金属丝半径；r′为沿金属丝半径方向距离。

2　定压下金属丝对推进剂传热过程仿真

大量的实验数据显示，嵌金属丝端燃装药，在出现压强峰后，最终压强会达到稳定，说明在其他参数如金属丝直径、装药物性都确定的情况下，对于定压强，沿金属丝装药燃速会达到一个稳定的数值，因为如果沿金属丝装药燃速没有达到稳定，则燃面和压强都不可能达到稳定平衡态。选择不同定压情况，进行金属丝对推进剂传热过程仿真，探究沿金属丝推进剂燃速、增速比与压强的关系，为内弹道计算和分析建立基础。

2.1仿真计算模型

仿真计算模型：金属丝采用银丝，直径为0.36 mm，熔点为1 235 K；装药初始温度298 K，初温燃速系数a=1.272 8×10-4，压强指数为n=0.35。火焰温度3 600 K，燃烧室压强取定值，气体常数R=287.43 J/(mol·K)。相关材料参数如表1。

由于银丝对燃面上装药温度的影响范围很小，径向影响范围在0.1 mm量级。为减小计算量，传热计算模型取金属丝附近0.3 mm半径范围的装药，并假设在该极小范围内燃面是平的，没有倾角。定压下传热计算与整个装药燃面的变化过程是无关的，只与金属丝附近装药推移有关系。因此，只需要建立耦合金属丝附近装药推移过程的传热分析，就可计算沿金属丝燃速，传热分析模型如图2所示。热传导分析可计算出温度分布，按照式(6)计算沿金属丝装药温度，并由此得到装药燃速。

表1　计算模型物性参数

图2　传热仿真模型示意图

2.2仿真方法

通过python编程对有限元软件进行二次开发，实现金属丝对装药的传热仿真计算，程序计算步骤如下：

(1)建立包含金属丝和装药的参数化传热计算模型。

(2)在Δt时间内进行传热仿真计算。

(3)程序提取结果文件中药柱燃面温度Ti，计算装药燃速系数和燃速ri，燃面推移距离Δd=riΔti。确定金属丝到达熔点的单元。

(4)通过脚本参数化方法改变计算模型，等同于去除沿金属丝装药燃烧的单元和金属丝熔断单元，并将温度场插值到新构型下，作为此构型下的初始温度场，更新传热边界。

(5)回到步骤(2)，不断地循环计算，直到第i+1步的沿银丝推进剂燃速ri+1与第i步燃速ri差值很小时，则认为燃速达到平衡，结束传热仿真计算。

2.3数值计算结果

通过热传导和金属丝处装药燃面推移仿真计算过程，得到了不同压强下沿金属丝装药温度和燃速随时间变化规律。

图3显示了沿金属丝装药温度-时间变化规律。装药初始温度为298 K，由于高温燃气提高了金属丝的温度，金属丝向附近装药进行加热，沿金属丝装药温度迅速上升。伴随着传热过程，装药燃烧，燃面不断的推移，且沿金属丝推进剂温度越高，燃速越大，推移越快。装药燃面推移速度的加快，会使得温度提升变缓，且当推移速度达到一定程度时，燃面上被加热的装药温度就无法继续上升，最终将达到稳定的平衡态。压强从1～5 MPa变化时，沿金属丝装药的平衡态温度从466 K降低到440 K。

沿金属丝装药燃速rw主要受装药温度和压强控制，装药温度越高，燃速系数越大，燃速也越大，由于温度随时间不断增加，所以燃速也随时间单调递增。但装药燃速和温度又相互影响，温度达到平衡时，燃速也同时达到平衡态，见图4中曲线。压强从1～5 MPa变化时，沿金属丝装药燃速从42.1 mm/s增大到52.9 mm/s。

图3　装药温度-时间关系曲线

图4　燃速rw-时间关系曲线

通过多个定压情况下的仿真分析，可计算出沿金属丝装药燃面上平衡态温度、燃速与压强的关系。图5、图6中曲线分别显示了沿金属丝装药燃面上平衡态温度随压强单调递减，沿金属丝推进剂燃速随压强单调递增。压强越高，推移越快，相当于缩短了金属丝对装药加热的时间，则其平衡态温度将有所降低。沿金属丝燃速由于受燃速系数(燃速系数减小)和压强(压强增大)两个因素影响，仍是增大的。

图5　平衡温度-压强曲线

图6　沿金属丝装药燃速-压强曲线

图7中，增速比随压强单调递减。增速比满足：

(7)

从式(7)可见，增速比等于燃速系数的比值，基础燃速的燃速系数是不变的，而沿金属丝装药温度随压强递减，则燃速系数比值随压强的增大减小。所以，增速比随压强增大也是减小的。

(8)

由此可得，在燃面变化过程中，金属丝附近锥角与增速比成反比；压强增大时，增速比减小，金属丝处燃面锥角变平缓；压强降低时，增速比增大，金属丝处燃面锥角变陡峭。该结论在下章中解释燃面变化过程中，将起到至关重要的作用。

图7　增速比-压强曲线

3　嵌金属丝端燃装药内弹道计算

嵌金属丝端燃装药在燃面推移过程中会出现初始压强峰，试验曲线如图8所示。

图8　试验压强-时间曲线

试验显示，压强峰是在点火后数秒时才出现的，且嵌金属丝装药燃面为锥形，流量较小，压强峰不可能由侵蚀燃烧造成。同时，恒定增速比模型(保持燃面锥角不变)数值仿真结果显示，压强随时间增大到最大值后平衡，没有出现达到峰值后回落的现象，因此，增速比为恒定数值的假定也是不能成立的。

采用python语言编程耦合传热过程和内弹道计算，通过传热计算得到沿金属丝装药燃速、增速比以及金属丝处燃面锥角，传入到内弹道计算程序代码中，进行燃面推移，提取燃面推移后燃面面积，并得到相应的压强，然后回代到传热计算程序中，再次计算沿金属丝燃速，通过增量耦合计算能够得到压强-时间曲线。仿真计算结果与试验结果吻合较好。

3.1计算模型

内弹道仿真计算模型：推进剂半径为60 mm，初温298 K，初温燃速系数a=1.272 8×10-4，压强指数为n，中心处嵌一根0.36 mm银丝，燃气温度3 600 K，喉部半径12 mm，气体常数R=287.43 J/(mol·K)。

模型计算重点考虑金属丝对内弹道燃面变化的影响，为简化模型，计算从端燃(无金属丝时)平衡态开始。

通过python控制有限元软件cae中的点、线的生成和移动，可实现燃面的自动推移过程，燃面推移方法如图9所示。由于金属丝处锥角会不断地变化，因此只能用折线近似模拟连续的燃面。金属丝处锥角决定了燃面的变化趋势，离开金属丝的燃面按平行推移。计算过程中在增量足够小、点数很多的情况下，能够达到较高的精度。

图9　燃面推移方法示意图

3.2内弹道仿真计算结果

采用增量计算方法，每两个燃面之间时间间隔为0.03 s，由于计算燃面过于密集，为便于显示，每隔10步输出一次结果，燃面变化过程如图10所示。

图10　燃面变化图

由燃面变化图可看出，燃面从开始的小锥面逐渐扩展，面积不断增大，到边界时，形状为一个外凸的锥面，随着燃面的变化，最终燃面演化为母线为直线的圆锥面，并达到平衡态。平衡态增速比为2.03，锥面角度约为60 °。

燃面面积随时间变化规律见图11，由图11中曲线可看到，燃面面积在0.4 s之前基本没有变化，仅有微弱的增加，随后出现急剧上升段，且在4.92 s达到峰值后回落，在7 s左右出现一个极小值，之后微弱上升，并达到平衡。

压强随燃面变化，与燃面具有相同的规律，见图12。4.92 s峰值压强为4.01 MPa，平衡态压强为3.34 MPa，若不嵌金属丝，平衡压强只有1.12 MPa，嵌金属丝后整个压强增加了3倍。对比图8试验曲线和图12中仿真得到的压强-时间曲线，可看到数值计算结果与试验结果压强变化规律具有很好的一致性。

图11　燃面面积-时间曲线

图12　压强-时间曲线

金属丝处燃面锥角是由沿金属丝燃速与基础燃速的比值即增速比决定的。增速比-时间关系曲线见图13。增速比在初期很短的时间内增大到最大值，而后下降，在4.92 s时达到一个极小值，最终达到平衡态。说明沿金属丝锥角开始由平缓急剧变尖，但随着压强增大，角度开始变大、变平，最后形成一个固定角度。

为说明燃面变化规律，选取了几个具有代表意义的燃面，作出了燃面变化示意图，如图14所示。

图13　增速比-时间关系曲线

图14　燃面变化示意图

初期锥面形成过程：金属丝对其附近装药进行加热，使得沿金属丝装药燃速较快，在金属丝处将形成锥角，随着平行推移过程，锥面将不断的扩大，燃面面积增大，压强也随之提高。压强的增大，导致增速比rw/r降低，即金属丝处锥面角度会变平缓。整个锥面母线是外凸的， 见图14中初期燃面。

最大燃面形成过程：由于锥面向边界推移和扩展，所以整个燃面面积不断增大，最后锥面扩展到边界附近时，燃面面积和压强都将达到最大值，见图14中最大燃面。此时，燃面特点是金属丝处燃面平缓，而边界处燃面陡峭。

燃面减小过程：锥面到达边界后，由于燃面是外凸的，随着燃面推移，金属丝附近区域平缓的燃面会扩大，而边界处斜度较大的燃面会在平行推移过程中迅速缩小并消失， 造成整个燃面面积迅速下降，压强也随之降低，最后燃面和压强都达到极小值，见图14中减小燃面。

燃面达到平衡过程：当压强达到极小值，增速比增大，金属丝处燃面锥角变陡，燃面面积开始缓慢增加，最终燃面在不断调整过程中逐渐变直，形成母线为直线的锥面，见图14中平衡燃面。随后的推移过程中，燃面面积不变。此时，燃速、增速比、压强均达到稳定平衡态。

4　结论

(1)在金属丝材质，直径和装药性质，燃烧室总温等参数确定的情况下，沿金属丝装药传热计算显示：压强增大时，沿金属丝装药燃速也随之增大，金属丝附近被加热的装药温度相对降低，增速比减小，金属丝处燃面锥角变平缓。

(2)内弹道仿真计算结果显示，嵌金属丝端燃装药初始压强峰是由于银丝对装药增速效果随燃烧室压强发生变化，即增速比和金属丝处燃面锥角变化引起的，燃面最大时，锥面是外凸的，燃面平衡时，压强稳定，燃面是一个母线为直线的锥面。

(3)本文所建立的嵌银丝端燃装药耦合传热过程的内弹道计算方法，成功地模拟了嵌金属丝端燃装药点火起始阶段的压强峰形成原因，为端燃装药点火过程优化提供了较准确的模型基础和计算方法。

[1]张炜,朱慧,夏智勋,等.嵌银丝贫氧推进剂燃烧特性[J].国防科技大学学报,2001,23(6):17-19.

[2]King M.Analytical modeling of effects of wires on solid motor ballistics[R].AIAA 89-2784.

[3]Gossant B,Godfroy F,Robert P H.Theoretical calculus of burning rate ratio in grains with embedded metal wires[R].AIAA 88-3255.

[4]张有为,王晓宏.嵌金属丝端燃药柱固体发动机水下点火数值模拟研究[J].固体火箭技术,2007,30(1):30-38.

[5]张有为,王晓宏.利用嵌金属丝药柱调节固体火箭发动机推力的计算研究[J].兵工学报,2007,28(10):1218-1222.

[6]张有为,王晓宏,杨举贤.嵌金属丝端燃药柱燃烧过程的数值研究[J].推进技术,2007,28(1):4-8.

[7]肖志斌,王家鑫,王继.嵌金属丝端燃装药方案绝热层设计方法[J].上海航天,2010(2):61-64.

[8]于宁,杨涛,刘巍,等.嵌银丝端燃药柱燃气发生器理论性能分析[J].弹箭与制导学报,2009,29(1):180-183.

[9]肖志平,张泽远,邢国强.星型药柱固化降温瞬态热力耦合分析[J].航空兵器,2014(1):40-43.

[10]王勖成.有限单元法[M].北京:清华大学出版社,2003.

[11]杨世铭,陶文铨.传热学(4版)[M].北京:高等教育出版社,2006.

[12]董师颜,张兆良.固体火箭发动机原理[M].北京:北京理工大学出版社,1996.

(编辑：崔贤彬)

Numerical study on the interior ballistics of grain embedded with silver-wire coupled convective heat transfer process

XIAO Zhi-ping，HAN Bo，ZHANG Liang，ZHI Shi-jun，XING Guo-qiang

(Air-to-Air Missile Research Institute，Luoyang471000，China)

The study focuses on the mechanism of the appearance of initial pressure peak when using the embedded wire on the end burning grain.The convective heat transfer process was modeled and the computational model for the velocity of burning was created.The heat transfer process and the motion of the burning grain were simulated via the secondary development of the finite element software.The result indicates that the temperature,the burning velocity and augmentation ratio of the grain along the wire is strongly affected by the pressure in the combustion chamber.Based on this finding,the interior ballistics simulation of the grain with only one embedded wire was simulated ,thereby the pressure ,the burning area and the burning velocity along the wire versus time were obtained.The good agreement with the experimental data indicates that the computational model can be used properly to simulate the burning process of the grain with wires. In addition,the appearance of the initial pressure peak during the burning process was explained theoretically.

embedded wire grain；secondary development；interior ballistics；heat transfer；burning rate

2015-04-28；

2015-06-15。

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

肖志平(1986—)，男，工程师，研究方向为发动机设计。E-mail：zhiping120@126.com

V435

A

1006-2793(2016)02-0188-06

10.7673/j.issn.1006-2793.2016.02.007