龙格库塔法计算固体火箭发动机内弹道
固体火箭发动机内弹道学 方丁酉
固体火箭发动机内弹道学方丁酉固体火箭发动机内弹道学方丁酉近年来，随着航天事业的蓬勃发展，固体火箭发动机内弹道学成为越来越受关注的领域。
方丁酉，中国工程院院士、火箭技术专家，曾在这一领域做出了卓越的贡献。
本文将从简单到复杂、由浅入深地探讨固体火箭发动机内弹道学，以帮助读者全面、深入地理解这一话题。
一、固体火箭发动机内弹道学的概念固体火箭发动机内弹道学是研究固体火箭发动机内燃烧过程、燃料燃烧特性、燃烧产物排放、工作环境等方面的学科。
方丁酉指出，固体火箭发动机内弹道学是固体火箭技术中的重要基础理论，对于提高固体火箭发动机的性能、可靠性和运载能力具有重要意义。
二、固体火箭发动机内弹道学的理论基础在固体火箭发动机内弹道学的研究中，燃烧动力学是重要的理论基础之一。
燃烧动力学研究了燃料在燃烧过程中的变化规律，以及燃烧反应对火箭发动机内部气体流动和压力变化的影响。
燃烧产物排放和燃烧室内部气体动力学也是固体火箭发动机内弹道学的重要内容之一。
研究这些理论基础可以帮助我们更好地理解固体火箭发动机内部的工作原理和特性。
三、固体火箭发动机内弹道学的关键技术在固体火箭发动机内弹道学研究中，燃烧稳定性和效率是两个关键技术。
燃烧稳定性是指在固体火箭发动机工作过程中保持燃烧的稳定性，避免出现燃烧不均匀或燃烧失稳等问题。
而燃烧效率则是指在燃烧过程中尽可能提高燃料的利用率，减少燃料的浪费。
方丁酉在固体火箭发动机内弹道学的研究中，提出了一系列有效的技术方案，使固体火箭发动机在燃烧稳定性和效率方面取得了显著的进展。
四、固体火箭发动机内弹道学的应用前景固体火箭发动机内弹道学的研究成果已经在我国的航天事业中得到了广泛的应用。
在长征系列火箭、嫦娥探月工程等多个航天工程中，固体火箭发动机内弹道学的研究成果为提高火箭的性能和可靠性做出了重要贡献。
未来，随着我国航天事业的不断发展，固体火箭发动机内弹道学将继续发挥重要作用，推动我国航天事业迈向新的高度。
龙格库塔法解轨道参数
龙格库塔法解轨道参数引言龙格库塔法（Runge-Kutta method）是一种常用的数值计算方法，用于求解常微分方程的数值解。
在天体力学中，我们经常需要通过数值方法来计算天体的轨道参数，如轨道椭圆的长短轴、离心率、倾角等。
本文将介绍龙格库塔法在解轨道参数中的应用，并详细探讨该方法的原理和实现过程。
基本原理龙格库塔法是一种迭代求解的方法，在每个时间步长内利用当前的状态来估计下一个状态。
具体而言，龙格库塔法将微分方程的求解问题转化为一个迭代的求解问题，通过逐步迭代来逼近精确解。
在解轨道参数的问题中，我们通常需要根据已知的初始条件以及天体的质量和力学模型来求解天体的轨道参数。
常用的力学模型有开普勒模型和牛顿模型。
龙格库塔法可以根据力学模型的不同进行相应的求解。
开普勒模型下的轨道参数求解步骤一：确定初始条件在使用龙格库塔法求解轨道参数之前，我们需要确定一些初始条件。
这些初始条件包括天体的质量、位置和速度。
步骤二：选择时间步长在求解过程中，我们需要选择一个合适的时间步长。
时间步长越小，计算的精度会越高，但计算的时间会增加。
步骤三：迭代求解利用龙格库塔法进行迭代求解的具体步骤如下：1.根据当前时刻的位置和速度，计算天体在该时刻的加速度。
2.根据当前时刻的位置、速度和加速度，计算下一个时刻的位置和速度。
3.更新当前时刻的位置和速度为新的位置和速度。
4.重复上述步骤，直到达到指定的终止条件。
步骤四：计算轨道参数通过迭代求解，我们可以得到天体在不同时刻的位置和速度。
根据这些位置和速度，我们可以计算出轨道参数，如离心率、倾角、长短轴等。
常用的轨道参数计算公式如下：1.离心率：e=√1+2El2μ(GM⊕)2)2.倾角：i=arccos(ℎzℎ3.长轴：a=−μT22E4.短轴：b=a√1−e2其中，E表示能量，l表示轨道角动量，μ表示标准引力参数，G表示引力常数，M⊕表示地球的质量，ℎ和ℎz分别表示轨道角动量和轨道角动量在z轴上的分量。
论业余模型火箭发动机设计方法5.31(1)
CH4
2.91E-07
2.77E-07
4.03E-07
CO
1.91E-01
1.85E-01
1.62E-01
CO2
1.54E-01
1.57E-01
1.75E-01
H
4.34E-06
2.38E-06
燃料名称 氧化剂 粘结剂
理论最大比冲
KNSU KNO3
蔗糖
137
KNSB KNO3 山梨醇
164
KNDX KNO3 葡萄糖
160
燃速数据如下
表2
图3
对于上述的燃料，可通过 Cprepop 软件进行相关参数的求解，下图为 Cprepop 软件界面
4
中青为科技小组
论业余模型火箭发动机设计方法
图4
下面以 KNSU 为例，在软件界面中输入相应的配比，这里硝酸钾对蔗糖之比为取 13:7。 可以计算化学平衡流与化学冻结流下的燃料产物组分与相关参数，如燃烧温度 Tf ，燃
固体火箭发动机的特点是通过消耗自身携带的燃料，向后排出高速运动的工质推动 载荷进行运动。一般而言，模型火箭要求其发动机制造价格低廉，结构简单，可快速重 新装填入模型火箭中进行发射作业。故现对一般的固体火箭发动机结构给予适当简化。 模型火箭发动机典型工作时间通常在 1.0s 左右，所以对于火箭发动机内壁的热防护可以 基本不必考虑，隔热层并没有安装的必要。对于市面上大多数模型固体火箭发动机而言， 通常都采用简单的圆孔尾喷管代替一般固体火箭发动机的拉法尔喷管；采用简单的深/ 浅内孔燃烧的管状装药以方便批量生产。可以说，模型火箭牺牲了一部分性能满足了上 述的价格低廉，结构简单的特点。图 1 所示为一般模型火箭发动机结构：
2 模型固体火箭发动机一般设计流程
内弹道 龙格库塔 计算 matlab
内弹道是指射程较短的导弹或火箭弹在飞行过程中受到大气阻力和重力等作用的飞行轨迹。
内弹道理论研究的是导弹或火箭弹在发射后到离开大气层再进入大气层末时的飞行过程。
内弹道包括导弹或火箭弹在发射后的加速、稳定、制导、飞行以及飞行过程中的动力学性能仿真等诸多内容。
内弹道有着复杂的飞行特性和动力学方程，在实际工程中需要进行准确的计算和仿真。
内弹道的计算中，龙格库塔（Runge-Kutta）法是一种常用的数值积分方法，在求解微分方程等领域有着广泛的应用。
龙格库塔法是由数学家奥特翁格（C. W. Runge）和马丁庫塔（M. W. J. Kutta）于1900年提出的，用于求解常微分方程初值问题，其优点是精度较高，适用范围广。
在内弹道计算中，可以利用龙格库塔法对导弹或火箭弹的飞行轨迹进行数值模拟和计算，得到较为准确的飞行轨迹数据。
在实际工程中，为了方便进行内弹道的计算，可以使用Matlab等数学建模和仿真软件。
Matlab是一种常用的科学计算软件，具有强大的数值计算和仿真功能，可以用于内弹道计算中的龙格库塔法数值模拟。
在Matlab中，可以编写相应的程序，利用龙格库塔法对导弹或火箭弹的飞行过程进行仿真和计算，得到准确的飞行轨迹和动力学性能数据。
内弹道计算是导弹或火箭弹研究设计中的重要内容，龙格库塔法是一种常用的数值积分方法，Matlab是一种常用的科学计算软件，它们的应用能够有效地进行内弹道的计算和仿真，为导弹或火箭弹的研制提供重要的技术支持。
随着技术的不断发展，内弹道计算已经成为导弹或火箭弹研究设计中不可或缺的一部分。
在内弹道计算中，龙格库塔法是一种常用的数值积分方法，可以对导弹或火箭弹的飞行轨迹进行数值模拟和计算，提供准确的飞行轨迹数据。
而Matlab作为一种强大的科学计算软件，对于内弹道的计算和仿真也有着重要的应用价值。
在实际工程中，使用Matlab编写程序，利用龙格库塔法对导弹或火箭弹的飞行轨迹进行数值模拟和计算，将为导弹或火箭弹的研制提供重要的技术支持。
固体火箭发动机燃烧室与内弹道计算
t
x
x
.
(V
V x
dx)dt
pr
Ab x
dxdt
m
整理后有：
t
(
A)
x
(
AV
)
pr
Ab x
燃烧面积
.. m+dm
. dm
dm
p
r
Ab x
dx
（2）动量方程
根据动量守恒：在dt时间内，微元体中动量的变化量=燃气 流动迁移在微元体内产生的动量增量+作用于两截面上的总压力 冲量+作用于微元体侧表面上压力冲量在x轴上的分量，即
x
(Ec
V2 2
)dx
dt
x
(
pAV
)dxdt
pr
Ab x
dxdt
Ip
整理后有：
t
A(Ec
V2 2
)
x
AV
(Ic
V2 2
)
prI
p
Ab x
Ec : 单位质量燃烧产物的内能(化学内能+热内能) Ic : 单位质量燃烧产物的总焓
1. 流动特点
燃烧室中的燃气流动具有以下特点： 燃烧室中推进剂装药不断燃烧，产生新的燃气加入主流，因此燃 气在燃烧室中的流动是一个有质量加入的流动过程； 在推进剂装药燃烧过程中，燃气的通道面积不断增大，故燃气的 流动参数是位置与时间的函数，因而是非定常流动； 燃气流过形状复杂的通道或流线变化剧烈时，将产生一定的能量 损失； 燃气的流动特性与推进剂的燃烧特性密切相关。例如燃气的压强 和流速影响推进剂的燃速，而燃速又反过来影响燃气压强和流速； 由此可见，在发动机工作过程中，燃气在燃烧室中的流动中极其 复杂的。为了研究方便，特作如下假设： 推进剂燃烧在装药燃烧表面上瞬时完成，燃气的化学组分和热力 性质均匀一致，并沿燃烧表面的外法线方向注入通道。 燃气为完全气体，服从完全气体状态方程，而且比热不变。 燃气在装药通道中作一维流动。 燃气与外界无热和功的交换。
【论文】龙格库塔法的推广及其应用
摘要在日常生活和工作中，我们经常遇到很多问题。
例如机械、电讯、核能、火箭、人造卫星、生物、医学及若干社会学科（如人口理论、经济预测等）的各个领域的问题，尤其是弹道轨道的定位、大型机械振动的分析、自动控制的设计、气象数值预报、按龄人口增长宏观预测问题等等。
这些问题都归结于求解常微分方程。
常微分方程的研究与应用已经深入到自然科学和社会科学的众多领域，并且成功地揭示了许多自然和社会现象的内在规律。
数值方法是求解常微分方程的一种有效的方法。
数值方法有很多种，龙格库塔法是解决初值问题的一种有效的数值方法。
我们可以通过选取不同参数从而得到不同的龙格库塔格式。
对于一类特定的初值问题，我们可以用选取恰当参数的龙格库塔格式求解这一类特定的初值问题。
每一种方法都有其自身的优点和缺点。
通过与解决这一类特定的初值问题的欧拉法的分析比较，我们可得到这两种方法的优点和缺点。
关键词：常微分方程，龙格库塔法，初值问题，导弹动力学弹道方程组IIIAbstractIn daily life and work, we often encounter many problems. Such as machinery, telecommunications, nuclear energy, rockets, satellites, biotechnology, medicine, and number of social disciplines (such as population theory, economic forecasts, etc.) in all areas, Especially ballistic orbit location, large machinery vibration analysis, control design, numerical weather forecasting, macroeconomic forecasts by age problems like population growth. These problems are due to solution of ordinary differential equations. Research and application of ordinary differential equations have been deep into many natural and social sciences fields, and successfully reveals the many natural and social phenomena of the inherent laws. Numerical solution of ordinary differential equations is an effective method. There are many numerical methods, Runge-Kutta method to solve initial value problems is an effective numerical method. We can choose different parameters resulting in different Runge-Kutta format. For a class-specific initial value problem, we can select the appropriate parameters of Runge-Kutta Schemes for this type of specific initial value. Each method has its own advantages and disadvantages. With the initial value to solve this particular type of problem analysis and comparison of the Euler method, we can obtain the advantages and disadvantages of both methods.Keywords：Ordinary Differential Equations, Runge-Kutta method, initial value problem,Ballistic missile dynamics equationsIV目录第一章前言1.1 常微分方程的应用 (1)1.2 数值方法的基本思想和内容 (1)1.3 龙格库塔法研究导弹轨迹的目的和意义 (3)第二章研究内容与方法2.1 选参构建龙格库塔格式 (4)2.2 误差分析 (7)2.3 简化导弹动力学方程组 (11)2.4 欧拉法求解导弹动力学方程组数值 (12)2.5 龙格库塔法求解导弹动力学方程组数值 (13)第三章研究结果与讨论3.1 欧拉法求解结果 (17)3.2 龙格库塔法求解结果 (17)3.3 两种方法求解结果的讨论 (17)第四章结论与展望4.1 结论 (19)4.2 改进与展望 (19)参考文献 (20)致谢 (21)I第一章前言1.1 常微分方程的应用求解常微分方程是解决现实生活中很多问题的有效的方法。
固体火箭发动机内弹道性能的仿真研究_刘宝华
2009年10月第10期电 子 测 试ELEC TRONIC T ES TOct.2009No.10固体火箭发动机内弹道性能的仿真研究刘宝华1,杨志菊2(1辽宁葫芦岛市海军飞行学院教研部 葫芦岛 125000,2辽宁葫芦岛92941部队 葫芦岛 125001) 摘 要:研究发动机燃烧室内压强随时间变化的规律,是固体火箭发动机工作过程分析的主要方面,一般采用实验法可以获得直观可靠的数据,但实验法耗资大、周期长,不易操作。
本文根据零维内弹道数学模型,运用龙格库塔法及M A T L A B语言对某型固体导弹发动机内弹道工作过程进行数值仿真,画出燃烧室内压强随时间变化的曲线,并进一步分析得出影响发动机内弹道性能的因素,仿真结果与发动机燃烧室内工作情况相符,影响因素分析与实验相一致,为固体火箭发动机内弹道性能研究提供了参考。
关键词:发动机;内弹道;数值仿真;燃烧室中图分类号:T N957 文献标识码:BSimulate research of solide rocket eng ine inside trajecto ry characteristicLiu Baohua1,Y ang Zhiju2(1T eaching and scientific resea rch ministry,N aval Flying A czdemy,Huludao125001,China;292941U nit,PL A,H uludao125001,China)A bstract:Researching the rule of the pressure with the time changing o f the eng ine's firebox is the m ost aspect of analy zing the process o f the solid missile engine.Co mmo nly,experimenta-tion can achieve the intuitionistic and credible data,but this method consume m uch tim e and money,and not easy to manipulate.According to zero-dimensional inside trajecto ry numerical simulate model,to o ur co untry some ty pe solid missile engine,applacating Lo nger-Kuta me thod and the lang uage of m atlab launches the inside trajecto ry numerical simulate of the w o rking pro cess,and draw ed the burner pressure-time curve,passing further analy sis o uted the influence facto rs of inside trajecto ry characteristic.The result of simulatio n is consistent to the fact in en-gine's firebox and analysis of the influence facto rs is consistent to experim entation.It offer s a new idea to researching the capability o f the ballistic trajecto ry.Keywords:engine;inside trajectory;numerical sim ulate;firebox0 引 言固体火箭发动机内弹道学的核心是研究发动机燃烧室内压强随时间变化的规律,是固体火箭发动机工作过程分析的主要方面。
第四讲龙格-库塔方法
龙格-库塔方法3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的一般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值方法，若用显式单步法(3.2.2)当，即数值求积用左矩形公式，它就是Euler法(3.1.2)，方法只有一阶精度，若取(3.2.3)就是改进Euler法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数f的值，但方法是二阶精度的.若要得到更高阶的公式，则求积分时必须用更多的f值，根据数值积分公式，可将(3.2.1)右端积分表示为注意，右端f中还不能直接得到，需要像改进Euler法(3.1.11)一样，用前面已算得的f值表示为(3.2.3)，一般情况可将(3.2.2)的 表示为(3.2.4)其中这里均为待定常数，公式(3.2.2)，(3.2.4)称为r级的显式Runge-Kutta法，简称R-K方法.它每步计算r个f值(即)，而k由前面(i-1)个已算出的表示，故公式是显式的.例i如当r=2时，公式可表示为(3.2.5) 其中.改进Euler 法(3.1.11)就是一个二级显式R-K 方法.参数取不同的值，可得到不同公式.3.2.2 二、三级显式R-K 方法对r=2的显式R-K 方法(3.2.5)，要求选择参数，使公式的精度阶p 尽量高，由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6) 令，对(3.2.6)式在处按Taylor 公式展开，由于将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须(3.2.7)即由此三式求的解不唯一.因r=2，由（3.2.5）式可知，于是有解(3.2.8)它表明使(3.2.5)具有二阶的方法很多，只要都可得到二阶精度R-K方法.若取，则，则得改进Euler法(3.1.11)，若取，则得，此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)及中点公式(3.2.9)是两个常用的二级R-K方法，注意二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个方程,加上(3.2.7)的三个方程，共计六个方程求4个待定参数，验证得出是无解的.当然r=2,p=2的R-K方法(3.2.5)当取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.对r=3的情形，要计算三个k值，即其中将按二元函数在处按Taylor公式展开，然后代入局部截断误差表达式，可得可得三阶方法，其系数共有8个，所应满足的方程为(3.2.10)这是8个未知数6个方程的方程组，解也是不唯一的，通常.一种常见的三级三阶R-K方法是下面的三级Kutta方法：(3.2.11)附：R-K 的三级Kutta 方法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满足c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K 方法及步长的自动选择利用二元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K 方法，其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =，2221(,(,))n n n n k f x a h y b hf x y =++，而33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定，Taylor 展开分析局部截断误差，使得精度达到四阶，即误差为5()O h 。
双脉冲火箭发动机设计文献综述
---------------------------------------------------------------范文最新推荐------------------------------------------------------ 双脉冲火箭发动机设计+文献综述摘要20世纪50年代，固体火箭在技术上取得突破，并广泛应用于弹道导弹和运载火箭。
20世纪60年代起，国外在可控推力固体火箭发动机的理论和实验研究方面作了大量的工作，探索出了很多技术途径和设计方案。
固体火箭发动机具有使用安全性好、可靠性高、储存性能好、密度比冲高及勤务处理方便等优点，使其在战略、战术导弹武器领域内成为主要的动力装臵。
脉冲固体火箭发动机是兴起于20世纪末期的新型固体火箭发动机，它的出现为推进系统的发展开辟了新的发展方向，在航天和军事领域都有着潜在的应用价值。
本论文主要研究双脉冲火箭发动机总体结构、装药结构设计、燃烧室设计、喷管设计、点火装臵设计以及内弹道计算。
通过对双脉冲固体火箭发动机的设计了解了双脉冲固体火箭发动机的结构以及关键技术，为研制高性能固体火箭发动机奠定坚实的基础。
60131 / 16关键词脉冲固体火箭发动机动力装臵总体结构装药结构毕业设计说明书（论文）外文摘要TitleDouble Pulse Solid Rocket Motor DesignAbstractIn the 1950's, the solid rocket breakthrough a lot in the technical field , and the solid rocket motor are widely used in ballistic missiles and launch vehicles. In the 1960's, the foreign countries made a great mass of works in theoretical and experimental studies of controllable thrust solid rocket motor to explore a lot of technical approaches and designs. Solid rocket motor with the advantage of security, reliability, and storage performance density than the ascribed and service processing advantages of convenience, it become a major power plant in the strategic field of tactical missiles weapons. Pulse solid rocket motor is a new type of solid rocket motor which risen in the late 20th century, its---------------------------------------------------------------范文最新推荐------------------------------------------------------propulsion system development has opened a new direction of development, have potential applications in the aerospace and military fields. This paper mainly studies the dual-pulse rocket engine general structure, propellant structure design, the combustion chamber design, nozzle design, the ignition system design and interior ballistic calculation. The purpose of designing dual-pulse solid rocket motor is to understand the structure and key technology of dual-pulse solid rocket motor, and lay a solid foundation for the development of high-performance solid rocket motor.4．7燃烧室的强度校核255发动机喷管设计265．1喷管的型面设计265．2喷管壁厚273 / 165．3喷管的热防护276点火装臵的设计296．1点火药的选择296．2点火药量的计算296．3点火药盒设计307内弹道的计算317．1内弹道计算的基本方程31 7．2四阶龙格-库塔法介绍31 7．3计算步骤327．4曲线绘制33结论35---------------------------------------------------------------范文最新推荐------------------------------------------------------ 致谢36参考文献371绪论1．1课题的研究背景20世纪50年代，固体火箭在技术上取得突破，并广泛应用于弹道导弹和运载火箭。
固体火箭发动机零维两相内弹道研究
固体火箭发动机零维两相内弹道研究陈军【摘要】为方便应用两相内弹道流动模型对火箭发动机进行性能预示以及提高性能预示精度,利用火箭喷管内的两相流动性能计算公式,建立了零维两相内弹道模型,包括零维两相内弹道微分方程和平衡压强公式,给出了模型中涉及到的两相特性参数的计算方法.相比于一维两相内弹道模型,该零维内弹道模型简单且满足必要的精度,适于工程应用的快速估算.利用该模型对某远程火箭发动机进行了内弹道计算与分析,计算结果与实验数据吻合良好,表明该两相内弹道模型可以有效地降低纯气相模型引起的理论与实际之间的模型偏差,有利于快速计算固体推进剂火箭发动机的两相内弹道参数以及提高预示精度.%To predict internal ballistic properties with two-phase flow for convenience and to improve prediction precision of internal ballistics in solid rocket motor(SRM),a model of zero-dimensional two-phase internal ballistics was built.The model includes internal ballistic differential equation and corresponding equilibrium pressure formula,and the computational methods of twophase property parameters in the model were pared with one-dimensional twophase internal ballistic model,the model is very simple and more suitable to quick computation in project applications,and more accurate.Based on the model,an internal ballistic computation and analysis for a long-range solid rocket engine were carried out.The computed results accord well with experimental data.The zero-dimensional two-phase internal ballistic model can availably reduce theoretical error caused by single-gas-phase model.By the model,two-phase internal ballisticproperties can be quickly computed,and the prediction precision of internal ballistics in SRM can be improved.【期刊名称】《弹道学报》【年(卷),期】2013(025)002【总页数】5页(P39-43)【关键词】内弹道;两相流;固体火箭发动机;固体推进剂【作者】陈军【作者单位】南京理工大学机械工程学院,南京210094【正文语种】中文【中图分类】V435随着高能推进剂在火箭发动机中的普遍应用，两相流动对发动机性能的影响越来越受到重视。
《火箭发动机》 7 内弹道
1 1 n
M
1 1 n
1 1 n 1 1 M km 1 n p RT0
当Peq较小时,如Peq=9.8MPa时,可取ε=0，则此时有： Peq M
1 1 n
2. 影响平衡压强的因素： 1 * 从平衡压强的表达式 Peq 1- C p aKN ( æ ) / 1-n 可以看到：当 处于理想状态（即流量修正系数ψ=1、热损失修正系数χ=1与ε=0）时， 影响平衡压强的因素主要有两方面：一方面是推进剂的特性，如密度、 特征速度、燃速特性与压强指数；另一方面是发动机几何设计参数，如 面喉比等。 (1)推进剂特性对平衡压强的影响 推进剂特性对平衡压强有着决定性的影响。由于推进剂的种类很多， 设计中可根据要求选择推进剂，确定相应的平衡压强。例如，就常用的 推进剂而言，其密度ρp约为1600-1800kg/m3，其特征速度C*约为12001600m/s，其燃速则差别更大，速燃推进剂可高达每秒几十毫米，缓燃 推进剂只有每秒几毫米，甚至每秒1毫米。所有这些数据都直接影响发 动机的平衡压强。
在发动机实际工作过程中，燃气密度ρ远小于装药密度ρp， 因此，可忽略填充量，故微分方程又可简化为：
Vg dp mb mt RT0 dt
式中 其中
mb p Ab r p Ab apn ()
为沿装药全长的平均侵蚀比。
pAt pAt C* RT0
根据质量守恒原理，燃烧室内燃气生成率 mb 与燃气通过喷管 排出的质量流率 mt 之差应等于燃烧室内燃气质量变化率，即：
dmr m b mt dt
p : 推进剂装药密度
Ab : 装药燃烧面积
m b p rAb 其中： mr Vg
高速旋转固体火箭发动机内流场数值计算
高速旋转固体火箭发动机内流场数值计算郝雯;封锋;罗盟;谢爱元【摘要】为计算高速旋转固体火箭发动机内流场特性,采用RNG k-ε湍流模型,对某端面和内孔同时燃烧管状装药固体火箭发动机燃烧室-喷管统一内流场进行了三维流场数值计算.采用用户自定义函数UDF编程给出了质量入口边界.计算结果表明,高速旋转固体火箭发动机内流场切向速度分布较复杂,在燃烧室前封头呈现出了准Rankie组合涡分布,在喷管内呈现明显的准强迫涡特性.燃气切向速度在燃烧室和喷管喉部分别达到峰值.燃烧室前封头处燃气轴向速度变化剧烈,出现了2个峰值,并在近轴线处出现了轴向速度的负值.高速旋转对前封头与喷喉结构强度与热防护带来严重影响.【期刊名称】《固体火箭技术》【年(卷),期】2013(036)004【总页数】5页(P481-485)【关键词】高速旋转;固体火箭发动机;组合涡;切向速度;轴向速度【作者】郝雯;封锋;罗盟;谢爱元【作者单位】南京理工大学机械工程学院,南京210094;南京理工大学机械工程学院,南京210094;南京理工大学机械工程学院,南京210094;南京理工大学机械工程学院,南京210094【正文语种】中文【中图分类】V4350 引言火箭增程及底排-火箭复合增程是当今炮射弹药增程的两种主要方式。
旋转稳定型弹药飞行时的转速可达到20 000 r/min或以上。
这种高速旋转的工作条件显著改变发动机的装药燃烧规律［1］及发动机内弹道参数［2］，主要表现在缩短固体火箭发动机的燃烧时间、增加推力与压强、降低比冲［3－5］。
此外，高速旋转所致的强旋流动现象对发动机结构的热防护带来负面效应［6］。
文献［7］中，采用某固体火箭发动机进行高速旋转实验研究，最后拟合出了动态燃速随装药燃面相对半径和转速变化的经验式。
本文参照此文献中拟合的燃速公式，在Fluent软件基础上，使用用户自定义函数UDF，给出质量入口边界，对端面和内孔同时燃烧的管状装药固体火箭发动机进行了三维内流场数值计算。
旋转固体火箭发动机一维内弹道计算
旋转固体火箭发动机一维内弹道计算
张为华;曹泰岳;郭印诚
【期刊名称】《推进技术》
【年(卷),期】1994()3
【摘要】给出了带＂锥－柱形＂装药的旋转固体火箭发动机的一维内弹道计算方法，其燃速计算采用稳态分析方法，喷管流量计算采用等熵旋转喷管流的简化方法。
所得结果和旋转试车台上实测的实验结果吻合较好。
计算和实验表明，＂锥－柱形＂装药呈现较大的加速度效应。
文中提供的计算方法具有较高预示精度，可推广到其它药四旋转发动机的内弹道计算，对旋转发动机设计具有指导作用。
【总页数】5页(P28-32)
【关键词】固体推进剂;火箭发动机;内弹道
【作者】张为华;曹泰岳;郭印诚
【作者单位】国防科技大学航天技术系
【正文语种】中文
【中图分类】V435.11
【相关文献】
1.固体火箭发动机燃烧室结构径向振动与内弹道性能的耦合计算 [J], 李亮;鲍福廷;孟伟;和兴锁
2.串联双燃速固体火箭发动机一维内弹道计算 [J], 鲍福廷;赵飞;李进贤
3.锥柱形装药旋转火箭发动机的零维内弹道计算 [J], 曹泰岳;吴军
4.固体火箭发动机点火过程内弹道计算 [J], 李海涛;鲍福廷
5.固体火箭发动机通用药柱的一维准定常内弹道计算 [J], 田维平;张扬中;顾雪林因版权原因，仅展示原文概要，查看原文内容请购买。
固体火箭发动机内弹道学 方丁酉
一、引言固体火箭发动机内弹道学方丁酉，是一门深奥而又具有重要意义的科学。
固体火箭发动机是航天领域中不可或缺的一部分，其内弹道学方丁酉的研究对于提高火箭发动机的性能和可靠性具有重要意义。
在本文中，将深入探讨固体火箭发动机内弹道学方丁酉的重要性、原理和应用。
二、固体火箭发动机内弹道学理论的重要性固体火箭发动机内弹道学方丁酉，是对固体火箭发动机内部燃烧过程的运动状态、燃烧特性、瞬态过程等进行研究的学科。
通过对固体火箭发动机内弹道学方丁酉的研究，可以深入了解燃烧室内各种参数的变化规律，确定燃烧室内流场的性质，以及改善发动机的设计和性能等方面提供重要参考。
固体火箭发动机内弹道学方丁酉的研究对提高火箭发动机的推进效率、降低燃料消耗、提高火箭发射成功率等具有重要意义。
通过对固体火箭发动机内弹道学方丁酉的研究，可以进一步提高火箭的响应速度和控制精度，提高火箭的战场生存能力和作战效能。
固体火箭发动机内弹道学方丁酉的研究内容主要包括燃烧室内燃气运动的多相流动、燃烧过程的热力学性质、燃烧产物的组成及其对流场和燃烧室壁面的影响等。
通过对这些内容的研究，可以从微观和宏观两个层面全面了解固体火箭发动机内弹道学方丁酉的运行机理，进而进行针对性的改进和应用。
固体火箭发动机内弹道学方丁酉的原理在航天领域中有着广泛的应用。
在燃烧室内燃气运动的多相流动方面，可以通过优化燃烧室的结构和燃烧室内部的燃烧过程，提高火箭的推进效率；在燃烧产物的组成及其对流场和燃烧室壁面的影响方面，可以减少火箭发动机的燃料消耗，提高火箭的响应速度和控制精度等。
四、固体火箭发动机内弹道学方丁酉的个人观点对于固体火箭发动机内弹道学方丁酉这一学科，我个人认为其研究的深入和广泛应用，对于航天领域的发展和火箭技术的提升有着极其重要的作用。
只有通过对固体火箭发动机内弹道学方丁酉的深入研究，才能够更好地理解火箭发动机内部燃烧过程的运动状态、燃烧特性、瞬态过程等，从而更好地改进火箭的设计和性能，提高火箭的推进效率、降低燃料消耗、提高火箭发射成功率等。
龙格-库塔(Runge-Kutta)方法
( )
dy dy y′ = = f, y′′ = f x + f y = f x + ffy = F dx dx F F dy F F y′′′ = + = +f x y dx x y F = f xx + f x f y + ffyx = f xx + f x f y + ffxy x F 2 2 f = f(fxy + f y f y + ffyy ) = ffxy + ffy + f f yy y
y ( x) = e
但
x2
∫
x
0
e dt
t2
∫
x
0
e dt 难以求积
t2
ODE数值解的基本思想和方法特点 数值解的基本思想和方法特点
1. 离散化 级数、 用Taylor级数、数值积分和差商逼近导数， 级数 数值积分和差商逼近导数， 将 ODE转化为离散的代数方程 称差分方程 。 转化为离散的代数方程(称差分方程 转化为离散的代数方程 称差分方程)。
(ha2 ) + (ha3 ) +
f 2 2! x
2 2 2
(ha f ) +
2
2
K3 f + ha3 f x + h (a3 b32 ) f + b32 K2 f y 2! 2! + h2a3 (a3 b32 ) f + b32 K2 f xy f xx + h (a3 b32 ) f + b32 K2
2
Euler法 后退 法 ym+1 = ym + hK2 + O(h2 ) K2 = f ( xm + h, ym+1 )
龙格库塔法计算固体火箭发动机内弹道
!计?算?星?型í装痢?药?的?几?何?尺?寸?!参?数簓符?号?说μ明÷!-----------------------------------------------------------------------------------!n表括?示?星?角?数簓,num表括?示?将?推?进?剂á沿?肉╝厚?方?向ò分?为a几?等台?分?,m表括?示?选?择?压1力 值μ大洙?小?!d表括?示?外猘径?,len表括?示?长¤度è,thet表括?示?星?边?夹D角?,epsilon表括?示?角?度è系μ数簓!r表括?示?过y度è圆2弧?半?径?,r1表括?示?星?角?圆2弧?半?径?,R0表括?示?通 ?用?气?体?常￡数簓!l表括?示?药?柱ù的?特?征÷长¤度è,y0表括?示?初?始?特?征÷参?数簓,y1表括?示?燃?尽?特?征÷参?数簓!I0表括?示?总哩?冲?,F表括?示?推?力 ,Poc表括?示?燃?烧?室酣?的?工¤作痢?压1力!Isp表括?示?比括?冲?,density_p表括?示?密ü度è,k表括?示?比括?热è?比括?rspeed表括?示?燃?速ù，?pn表括?示?压1力 指?数簓!mpeff表括?示?有瓺效§装痢?药?量?,Cf表括?示?推?力 系μ数簓,Ctz表括?示?特?征÷速ù度è,At表括?示?喉í部?面?积y,a表括?示?燃?速ù系μ数簓!S表括?示?平?均ù燃?烧?面?积y,e1表括?示?平?均ù肉╝厚?,epsilon1表括?示?减?面?比括?epsilon2表括?示?增?面?比括?!foresmax表括?示?前°段?最?大洙?相à对?周ü边?长¤,backsmax表括?示?后ó段?最?大洙?相à对?周ü边?长¤!smin表括?示?最?小?相à对?周ü边?长¤,thet1表括?示?周ü边?长¤取?得?最?小?值μ时骸?的?星?边?夹D角?!Ap表括?示?初?始?通 ?气?面?积y,J表括?示?初?始?通 ?气?参?量?,eta表括?示?装痢?填?系μ数簓,Af表括?示?剩骸?药?面?积y,etaf表括?示?剩骸?药?系μ数簓!de表括?示?每?一?份 肉╝厚?的?长¤度è,Sa表括?示?燃?烧?面?积y数簓组哩?Apa表括?示?通 ?气?面?积y数簓组哩?!------------------------------------------------------------------------------------program mainimplicit nonereal(kind=8),parameter::Pi=3.14integer::n,num,i,mreal(kind=8)::d,len,thet,epsilon,r,r1real(kind=8)::l,y0,y1real(kind=8)::I0,F,Poc,Pe,R0real(kind=8)::Isp,density_p,k,rspeed,Pnreal(kind=8)::mpeff,Cf,Ctz,At,areal(kind=8)::S,e1,epsilon1,epsilon2real(kind=8)::foresmax,backsmax,smin,thet1real(kind=8)::Ap,J,eta,Af,etafreal(kind=8)::error,dereal(kind=8),allocatable::Sa(:),Apa(:)!读á入?所ù需è要癮所ù用?参?数簓值μopen(3,file="design_parameter.dat")read(3,*)I0!总哩?冲?read(3,*)F!推?力read(3,*)Poc!燃?烧?室酣?压1力read(3,*)Isp!比括?冲?read(3,*)density_p!推?进?剂á的?密ü度èread(3,*)k!比括?热è?比括?read(3,*)a!燃?速ù系μ数簓read(3,*)Ctz!读á取?特?征÷速ù度èCtzread(3,*)Pn!压1力 指?数簓read(3,*)Pe!读á取?喷?管ü出?口ú的?压1力 Peread(3,*)R0!读á取?通 ?用?气?体?常￡数簓R0read(3,*)d,r!读á取?装痢?药?直±径?d和í过y度è圆2弧?半?径?rread(3,*)epsilon2!读á取?增?面?比括╡psilon2read(3,*)r1!读á取?星?角?圆2弧?半?径?r1read(3,*)thet1!试?取?周ü边?长¤取?得?最?小?值μ时骸?的?星?边?夹D角?thet1read(3,*)thet!试?取?初?始?时骸?的?星?边?夹D角?thetread(3,*)n!读á取?星?角?数簓n（辍?,4,5,6,7,8）?read(3,*)epsilon!试?取?角?度è系μ数簓epsilonread(3,*)num!读á取?等台?分?肉╝厚?的?等台?分?数簓numread(3,*)m!m=1表括?示?最?小?压1力 给?定¨,m=2表括?示?最?大洙?压1力 给?定¨,m=0表括?示?平?均ù压1力 给?定¨close(3)allocate(Sa(0:num),Apa(0:num))!根ù据Y规?定¨的?总哩?冲?计?算?有瓺效§装痢?药?量?mpeffmpeff=1.02*I0/Isp!计?算?推?力 系μ数簓和í喉í部?面?积yCf=sqrt(k)*(2/(k+1))**((k+1)/(k-1)/2)*sqrt(2*k*(1-(Pe/Poc)**((k-1)/k))/(k-1))At=F/Cf/(Poc*101325.0)!计?算?平?均ù肉╝厚?e1和í计?算?特?征÷长¤度è和í平?均ù燃?烧?面?积yS!若?给?定¨的?是?最?大洙?压1力 则ò需è根ù据Y增?面?比括?计?算?最?小?燃?面?再ù计?算?平?均ù燃?面?值μ!y1一?般?取?值μ在ú1附?近ü,(0.8-1.2)if(m==0)then!平?均ù压1力 给?定¨S=At*(Poc*101325.0)**(1-Pn)/(Ctz*density_p*a)e1=mpeff/(density_p*S)else if(m==1)then!最?小?压1力 给?定¨S=At*(Poc*101325.0)**(1-Pn)/(Ctz*density_p*a)S=(epsilon2+1.0)*S/2.0e1=mpeff/(density_p*S)else!最?大洙?压1力 给?定¨S=At*(Poc*101325.0)**(1-Pn)/(Ctz*density_p*a)S=(1.0/epsilon2+1.0)*S/2.0e1=mpeff/(density_p*S)end ifl=d/2-e1-ry0=r/ly1=(e1+r)/l!根ù据Yepsilon2选?择?角?度è系μ数簓epsilon!计?算?取?得?最?小?周ü长¤时骸?的?thet1do while(.true.)error=thet1/2+cotan(thet1/2)-Pi/n-Pi/2if(error>=0.0001)thenthet1=thet1+0.00001else if(error<=-0.0001)thenthet1=thet1-0.00001elseexitend ifend do!计?算?角?度è系μ数簓epsilon的?准?确ā?值μdo while(.true.)backsmax=2*n*((1-epsilon)*Pi/n+y1*(Pi/n+asin(sin(epsilon*Pi/n)/y1))) smin=2*n*(sin(epsilon*Pi/n)/sin(thet1/2)+(1-epsilon)*Pi/n)if(backSmax/Smin>epsilon2)thenepsilon=epsilon+0.001elseexitend ifend do!根ù据Y减?面?比括╡psilon1，?计?算?thetepsilon1=1/epsilon2foresmax=smin/epsilon1do while(.true.)error=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(r1+r)*(Pi/n+Pi/2-thet/2-cos(thet/2)/sin(thet/2))/l)-foresmaxif(error>=0.0001)thenthet=thet+0.00001else if(error<=-0.0001)thenthet=thet-0.00001elseexitend do!计?算?初?始?通 ?气?面?积y，?通 ?气?参?量?，?装痢?填?系μ数簓,药?柱ù长¤度èAp=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-&sin(epsilon*Pi/n)*cotan(thet/2)))+2*n*r*(sin(epsilon*Pi/n)/sin(thet/2)+&(1-epsilon)*Pi/n)/l+n*r**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2+&n*r1**2*(thet/2+cotan(thet/2)-Pi/2)/l**2)*l**2J=At/Apeta=4*(Pi*d**2/4-Ap)/(Pi*d**2)Af=(epsilon*Pi*(1+y1)**2-n*(sin(epsilon*Pi/n)*(sqrt(y1**2-sin(epsilon*Pi/n)**2)+&cos(epsilon*Pi/n)))-n*y1**2*(epsilon*Pi/n+asin(sin(epsilon*Pi/n)/y1)))*l**2etaf=4*Af/(pi*d**2)len=mpeff/density_p/(Pi*d**2/4-Ap-Af)!计?算?推?进?剂á的?燃?面?变?化ˉ规?律é并￠输?出?结á果?de=e1/numopen(10,file="export_burnS.dat")open(20,file="export_Ap.dat")do i=0,num,1if((i*de)<=r1)thenSa(i)=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(r1+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l-(r1-i*de)*Pi/n/l)*l*len Apa(i)=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-sin(epsilon*Pi/n)* cotan(thet/2)))+&2*n*(r+i*de)*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n)/l+&n*(r+i*de)**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2+n*(r1-i*de)**2*(thet/2+cotan(thet/2) -Pi/2)/l**2)*l**2else if((i*de)>r1.and.(i*de)<=(l*sin(epsilon*Pi/n)/cos(thet/2)-r))thenSa(i)=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(i*de+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l)*l*lenApa(i)=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-sin(epsilon*Pi/n)* cotan(thet/2)))+&2*n*(r+i*de)*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n)/l+&n*(r+i*de)**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2)*l**2else if((i*de)>(l*sin(epsilon*Pi/n)/cos(thet/2)-r).and.(i*de)<=e1)thenSa(i)=2*n*((1-epsilon)*Pi/n+(r+i*de)*(Pi/n+asin(l*sin(epsilon*Pi/n)/(i*de+r)))/l)*l *lenApa(i)=n*((1-epsilon)*Pi*(1+(r+i*de)/l)**2/n+sin(epsilon*Pi/n)*(sqrt((r+i*de)**2/l* *2-&sin(epsilon*Pi/n)**2)+cos(epsilon*Pi/n))+(r+i*de)**2*(epsilon*Pi/n+&asin(l*sin(epsilon*Pi/n)/(r+i*de)))/l**2)*l**2write(10,"(f10.5,2X,f15.5)")i*de,Sa(i)write(20,"(f10.5,2X,f15.5)")i*de,Apa(i)end doclose(10)close(20)!调獭?用?子哩?程ì序ò计?算?装痢?药?内ú弹獭?道台?曲ú线?callinternal_ballistics0(d,len,e1,n,thet,epsilon,r,r1,l,a,pn,Ctz,At,k,density_p,mpeff,Isp) !输?出?星?型í装痢?药?的?几?何?参?数簓open(30,file="export_star_geometry.dat")write(30,"(A6,f11.5,A2)")"mpeff=",mpeff,"Kg"write(30,"(A2,f8.5,A2)")"D=",d,"m"write(30,"(A2,f8.5,A2)")"L=",len,"m"write(30,"(A2,I2)")"n=",nwrite(30,"(A3,f8.5)")"θè=",thetwrite(30,"(A3,f7.5)")"ε?=",epsilonwrite(30,"(A2,f7.5,A2)")"r=",r,"m"write(30,"(A3,f7.5,A2)")"r1=",r1,"m"write(30,"(A3,f8.5,A2)")"e1=",e1,"m"write(30,"(A3,f10.5,A2)")"Ap=",Ap,"㎡O"write(30,"(A2,f10.5)")"J=",Jwrite(30,"(A3,f7.5)")"η?=",etawrite(30,"(A3,f10.5,A2)")"Af=",Af,"㎡O"write(30,"(A4,f10.5)")"η?f=",etafwrite(30,"(A2,f10.5)")"l=",lwrite(30,"(A3,f10.5)")"y0=",y0write(30,"(A3,f10.5)")"y1=",y1close(30)stopend!该?子哩?程ì序ò用?于 ?计?算?零?维?变?截?面?燃?烧?装痢?药?的?内ú弹獭?道台?!利?用?龙ⅷ?格?-库a塔t法ぁ?计?算?内ú弹獭?道台?曲ú线?!可é用?于 ?计?算?侵?蚀骸?燃?烧?效§应畖下?的?内ú弹獭?道台?曲ú线?!--------------------------------------------------------------------------!d外猘径?,len长¤度è,e1平?均ù肉╝厚?,n星?角?数簓,thet星?边?夹D角?,epsilon角?度è系μ数簓!r过y度è圆2弧?半?径?,r1星?角?圆2弧?半?径?和íl特?征÷尺?寸?,key表括?示?是?否?考?虑?侵?蚀骸?燃?烧?!density_p推?进?剂á的?密ü度è,k比括?热è?比括?ga系μ数簓,C特?征÷速ù度è,mpeff药?柱ù质ê量?!a速ù度è系μ数簓,pn压1强?指?数簓,At喉í部?面?积y,Pc燃?烧?室酣?的?设Θ?计?压1力 ,Isp表括?示?理え?论?比括?冲?!dt时骸?间?步?长¤,e燃?层?厚?度è,time时骸?间?,Sa同?一?时骸?刻ì轴á向ò的?各÷节ú点?的?燃?面?,Apa通 ?气?面?积y!Poc各÷节ú点?的?压1力 ,Pe喷?管ü出?口ú压1力 ,P_I压1力 冲?量?,ep侵?蚀骸?比括?!P_av平?均ù压1力 ,F_av平?均ù推?力 ,I0总哩?冲?,Im重?量?比括?冲?Iv体?积y比括?冲?,Cf 推?力 系μ数簓!--------------------------------------------------------------------------subroutineinternal_ballistics0(d,len,e1,n,thet,epsilon,r,r1,l,a,pn,Ctz,At,k,density_p,mpeff,Isp) implicit nonereal,parameter::Pi=3.14integer::n,keyreal(kind=8),intent(in)::d,len,e1,thet,epsilon,r,r1,lreal(kind=8),intent(in)::density_p,k,Ctz,mpeff,Ispreal(kind=8)::a,pn,At,dt,e,time,Sa,Apa,Poc,Pe,P_I,Freal(kind=8)::ep,P_av,F_av,I0,Im,Iv,Cf,gareal(kind=8)::f1,f2,f3,f4Pe=101325.0P_I=0.0!药?柱ù的?能ü量?特?性?参?量?读á入?open(3,file="inter_parameter.dat")read(3,*)keyread(3,*)dtread(3,*)Pocclose(3)!计?算?点?火e压1强?时骸?的?推?力 值μPoc=Poc*101325.0ga=sqrt(k)*(2/(k+1))**((k+1)/(k-1)/2)Cf=ga*sqrt(2*k*(1-(Pe/Poc)**((k-1)/k))/(k-1))F=Cf*Poc*Ate=0.0time=0.0open(10,file="export_Poc.dat")open(20,file="export_F.dat")write(10,"(f8.5,2X,f10.5)")time,Poc/101325.0write(20,"(f8.5,2X,f15.5)")time,Fdo while(.true.)!计?算?燃?面?和í通 ?气?面?积yif(e<=r1)thenSa=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(r1+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l-(r1-e)*Pi/n/l)*l*len Apa=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-&sin(epsilon*Pi/n)*cotan(thet/2)))+2*n*(r+e)*(sin(epsilon*Pi/n)/&sin(thet/2)+(1-epsilon)*Pi/n)/l+n*(r+e)**2*(Pi/n+Pi/2-thet/2-&cotan(thet/2))/l**2+n*(r1-e)**2*(thet/2+cotan(thet/2)-Pi/2)/l**2)*l**2else if(e>r1.and.e<=(l*sin(epsilon*Pi/n)/cos(thet/2)-r))thenSa=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(e+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l)*l*lenApa=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-&sin(epsilon*Pi/n)*cotan(thet/2)))+2*n*(r+e)*(sin(epsilon*Pi/n)/&sin(thet/2)+(1-epsilon)*Pi/n)/l+n*(r+e)**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2)*l**2 else if(e>(l*sin(epsilon*Pi/n)/cos(thet/2)-r).and.e<=e1)thenSa=2*n*((1-epsilon)*Pi/n+(r+e)*(Pi/n+asin(l*sin(epsilon*Pi/n)/(e+r)))/l)*l*lenApa=n*((1-epsilon)*Pi*(1+(r+e)/l)**2/n+sin(epsilon*Pi/n)*(sqrt((r+e)**2/l**2-& sin(epsilon*Pi/n)**2)+cos(epsilon*Pi/n))+(r+e)**2*&(epsilon*Pi/n+asin(l*sin(epsilon*Pi/n)/(r+e)))/l**2)*l**2elseSa=0.0end if!计?算?侵?蚀骸?比括?if(key==1)thenif(Sa/Apa<=72.9)thenep=1.0elseep=1.3128-1.3249e-2*Sa/Apa+1.5527e-4*(Sa/Apa)**2-4.3868e-7*(Sa/Apa)**3 end ifelseep=1.0end if!采é用?4阶×的?龙ⅷ?格?库a塔t法ぁ?计?算?f1=ga**2*Ctz**2*(density_p*Sa*ep*a*Poc**pn-Poc*At/Ctz)/(Apa*len)f2=ga**2*Ctz**2*(density_p*Sa*ep*a*(Poc+dt*f1/2)**pn-(Poc+dt*f1/2)*At/Ctz)/(Apa*len )f3=ga**2*Ctz**2*(density_p*Sa*ep*a*(Poc+dt*f2/2)**pn-(Poc+dt*f2/2)*At/Ctz)/(Apa*len )f4=ga**2*Ctz**2*(density_p*Sa*ep*a*(Poc+dt*f3)**pn-(Poc+dt*f3)*At/Ctz)/(Apa*len) Poc=Poc+dt*(f1+2.0*f2+2.0*f3+f4)/6.0Cf=ga*sqrt(2*k*(1-(Pe/Poc)**((k-1)/k))/(k-1))F=Cf*Poc*Attime=time+dte=e+dt*ep*a*Poc**pnif(e>e1.and.Poc<101325.0)exitwrite(10,"(f8.5,2X,f10.5)")time,Poc/101325.0write(20,"(f8.5,2X,f15.5)")time,FP_I=P_I+Poc*dtend doclose(10)close(20)!循-环·结á束?P_av=P_I/timeCf=ga*sqrt(2*k*(1-(Pe/P_av)**((k-1)/k))/(k-1))F_av=Cf*P_av*AtI0=F_av*timeIm=I0/mpeff/9.81Iv=I0/(Apa*len)!输?出?固ì体?火e箭y发ぁ?动ˉ机ú的?工¤作痢?参?数簓值μopen(30,file="export_propellant.dat")write(30,"('平?均ù压1力 p_av=',f9.5)")p_av/101325.0write(30,"('平?均ù推?力 F_av=',f10.5,1X,'KN')")F_av/1000.0 write(30,"('总哩?冲?I0=',f10.5,1X,'KN.s')")I0/1000.0write(30,"('重?量?比括?冲?Im=',f10.5,1X,'s')")Imwrite(30,"('体?积y比括?冲?Iv=',f10.5,1X,'KN.s/m3')")Iv/1000.0 write(30,"('工¤作痢?时骸?间?time=',f8.5,1X,'s')")timeclose(30)returnend subroutine。
《火箭发动机》 7 内弹道 共21页PPT资料
在上面的分析中，认为燃烧室是一个充满高压燃烧气体的容器， 不考虑燃气的流动和燃烧室内的压强分布，室内各点的压强都相等。 这样，整个燃烧室压强同时随时间变化，与该点的位置坐标x无关，这 就是所谓“零维”的压强变化。对于燃气流速很小的燃烧室来说，压 强计算可以看作是一个“零维”问题来处理。但是，对装填密度较大 的侧面燃烧装药，燃气在通道中的流动沿轴向产生很大的速度，因此， 压强沿轴向有显著的变化。这种情况下，必须考虑压强在燃烧室中的 分布，应作为“一维”问题来进行压强计算。
由发动机实验所测得的 燃烧室压强一时间曲线可见， 燃烧室压强的变化有三个阶 段，如右图所示：
1.发动机起动阶段(上升段) 这包括点火和压强建立过程。首先 依靠点火装置中点火药点燃并燃烧生成的高温气体充满燃烧室，一方 面使燃烧室压强上升到点火压强；另一方面加热推进剂表面，点燃主 装药，这就是点火过程。当主装药全面点燃后，燃气质量生成量迅速 增大，并在瞬时超过喷管的质量流量，使燃烧室的压强迅速增加，同 时又促使喷管流量的增加，不断地与燃气生成量趋于相对平衡。最后， 燃烧室压强达到其相对稳定值，这个相对稳定值的压强称为工作压强。 这个压强建立的过程即称为发动机启动阶段。对一般发动机来说，这 个过程在几十毫秒内完成。
第七章 固体火箭发动机的内弹道计算
一、内弹道计算的任务 二、燃烧室压强的变化 三、零维内弹道计算的微分方程 四、平衡压强及其影响因素 五、燃烧室压强—时间曲线的简化计算
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
!计?算?星?型³装痢?药?的?几?何?尺?寸?!参?数簓符?号?说μ明¶!-----------------------------------------------------------------------------------!n表括?示?星?角?数簓,num表括?示?将?推?进?剂®沿?肉╝厚?方?向´分?为a几?等台?分?,m表括?示?选?择?压1力 值μ大洙?小?!d表括?示?外猘径?,len表括?示?长¤度¯,thet表括?示?星?边?夹D角?,epsilon表括?示?角?度¯系μ数簓!r表括?示?过y度¯圆2弧?半?径?,r1表括?示?星?角?圆2弧?半?径?,R0表括?示?通 ?用?气?体?常￡数簓!l表括?示?药?柱·的?特?征¶长¤度¯,y0表括?示?初?始?特?征¶参?数簓,y1表括?示?燃?尽?特?征¶参?数簓!I0表括?示?总哩?冲?,F表括?示?推?力 ,Poc表括?示?燃?烧?室酣?的?工¤作痢?压1力!Isp表括?示?比括?冲?,density_p表括?示?密¹度¯,k表括?示?比括?热¯?比括?rspeed表括?示?燃?速·，?pn表括?示?压1力 指?数簓!mpeff表括?示?有瓺效§装痢?药?量?,Cf表括?示?推?力 系μ数簓,Ctz表括?示?特?征¶速·度¯,At表括?示?喉³部?面?积y,a表括?示?燃?速·系μ数簓!S表括?示?平?均·燃?烧?面?积y,e1表括?示?平?均·肉╝厚?,epsilon1表括?示?减?面?比括?epsilon2表括?示?增?面?比括?!foresmax表括?示?前©段?最?大洙?相对?周¹边?长¤,backsmax表括?示?后µ段?最?大洙?相对?周¹边?长¤!smin表括?示?最?小?相对?周¹边?长¤,thet1表括?示?周¹边?长¤取?得?最?小?值μ时骸?的?星?边?夹D角?!Ap表括?示?初?始?通 ?气?面?积y,J表括?示?初?始?通 ?气?参?量?,eta表括?示?装痢?填?系μ数簓,Af表括?示?剩骸?药?面?积y,etaf表括?示?剩骸?药?系μ数簓!de表括?示?每?一?份 肉╝厚?的?长¤度¯,Sa表括?示?燃?烧?面?积y数簓组哩?Apa表括?示?通 ?气?面?积y数簓组哩?!------------------------------------------------------------------------------------program mainimplicit nonereal(kind=8),parameter :: Pi=3.14integer :: n,num,i,mreal(kind=8) :: d,len,thet,epsilon,r,r1real(kind=8) :: l,y0,y1real(kind=8) :: I0,F,Poc,Pe,R0real(kind=8) :: Isp,density_p,k,rspeed,Pnreal(kind=8) :: mpeff,Cf,Ctz,At,areal(kind=8) :: S,e1,epsilon1,epsilon2real(kind=8) :: foresmax,backsmax,smin,thet1real(kind=8) :: Ap,J,eta,Af,etafreal(kind=8) :: error,dereal(kind=8),allocatable :: Sa(:),Apa(:)!读®入?所·需¯要癮所·用?参?数簓值μopen(3,file="design_parameter.dat")read(3,*) I0 !总哩?冲?read(3,*) F !推?力read(3,*) Poc !燃?烧?室酣?压1力read(3,*) Isp !比括?冲?read(3,*) density_p !推?进?剂®的?密¹度¯read(3,*) k !比括?热¯?比括?read(3,*) a !燃?速·系μ数簓read(3,*) Ctz !读®取?特?征¶速·度¯Ctzread(3,*) Pn !压1力 指?数簓read(3,*) Pe !读®取?喷?管¹出?口¸的?压1力 Peread(3,*) R0 !读®取?通 ?用?气?体?常￡数簓R0read(3,*) d,r !读®取?装痢?药?直ª径?d和³过y度¯圆2弧?半?径?rread(3,*) epsilon2 !读®取?增?面?比括╡psilon2read(3,*) r1 !读®取?星?角?圆2弧?半?径?r1read(3,*) thet1 !试?取?周¹边?长¤取?得?最?小?值μ时骸?的?星?边?夹D角?thet1read(3,*) thet !试?取?初?始?时骸?的?星?边?夹D角?thetread(3,*) n !读®取?星?角?数簓n（辍?,4,5,6,7,8）?read(3,*) epsilon !试?取?角?度¯系μ数簓epsilonread(3,*) num !读®取?等台?分?肉╝厚?的?等台?分?数簓numread(3,*) m !m=1表括?示?最?小?压1力 给?定¨,m=2表括?示?最?大洙?压1力 给?定¨,m=0表括?示?平?均·压1力 给?定¨close(3)allocate(Sa(0:num),Apa(0:num))!根·据Y规?定¨的?总哩?冲?计?算?有瓺效§装痢?药?量?mpeffmpeff=1.02*I0/Isp!计?算?推?力 系μ数簓和³喉³部?面?积yCf=sqrt(k)*(2/(k+1))**((k+1)/(k-1)/2)*sqrt(2*k*(1-(Pe/Poc)**((k-1)/k))/(k-1))At=F/Cf/(Poc*101325.0)!计?算?平?均·肉╝厚?e1和³计?算?特?征¶长¤度¯和³平?均·燃?烧?面?积yS!若?给?定¨的?是?最?大洙?压1力 则´需¯根·据Y增?面?比括?计?算?最?小?燃?面?再·计?算?平?均·燃?面?值μ!y1一?般?取?值μ在¸1附?近¹,(0.8-1.2)if(m==0) then!平?均·压1力 给?定¨S=At*(Poc*101325.0)**(1-Pn)/(Ctz*density_p*a)e1=mpeff/(density_p*S)else if(m==1) then!最?小?压1力 给?定¨S=At*(Poc*101325.0)**(1-Pn)/(Ctz*density_p*a)S=(epsilon2+1.0)*S/2.0e1=mpeff/(density_p*S)else!最?大洙?压1力 给?定¨S=At*(Poc*101325.0)**(1-Pn)/(Ctz*density_p*a)S=(1.0/epsilon2+1.0)*S/2.0e1=mpeff/(density_p*S)end ifl=d/2-e1-ry0=r/ly1=(e1+r)/l!根·据Yepsilon2选?择?角?度¯系μ数簓epsilon!计?算?取?得?最?小?周¹长¤时骸?的?thet1do while(.true.)error=thet1/2+cotan(thet1/2)-Pi/n-Pi/2if(error>=0.0001) thenthet1=thet1+0.00001else if(error<=-0.0001) thenthet1=thet1-0.00001elseexitend ifend do!计?算?角?度¯系μ数簓epsilon的?准?确º?值μdo while(.true.)backsmax=2*n*((1-epsilon)*Pi/n+y1*(Pi/n+asin(sin(epsilon*Pi/n)/y1))) smin=2*n*(sin(epsilon*Pi/n)/sin(thet1/2)+(1-epsilon)*Pi/n) if(backSmax/Smin>epsilon2) thenepsilon=epsilon+0.001elseexitend ifend do!根·据Y减?面?比括╡psilon1，?计?算?thetepsilon1=1/epsilon2foresmax=smin/epsilon1do while(.true.)error=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(r1+r)*(Pi/n+Pi/2-thet/2-cos(thet/2)/sin(thet/2))/l)-foresmax if(error>=0.0001) thenthet=thet+0.00001else if(error<=-0.0001) thenthet=thet-0.00001elseexitend do!计?算?初?始?通 ?气?面?积y，?通 ?气?参?量?，?装痢?填?系μ数簓,药?柱·长¤度¯Ap=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-&sin(epsilon*Pi/n)*cotan(thet/2)))+2*n*r*(sin(epsilon*Pi/n)/sin(thet/2)+&(1-epsilon)*Pi/n)/l+n*r**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2+&n*r1**2*(thet/2+cotan(thet/2)-Pi/2)/l**2)*l**2J=At/Apeta=4*(Pi*d**2/4-Ap)/(Pi*d**2)Af=(epsilon*Pi*(1+y1)**2-n*(sin(epsilon*Pi/n)*(sqrt(y1**2-sin(epsilon*Pi/n)**2)+& cos(epsilon*Pi/n)))-n*y1**2*(epsilon*Pi/n+asin(sin(epsilon*Pi/n)/y1)))*l**2etaf=4*Af/(pi*d**2)len=mpeff/density_p/(Pi*d**2/4-Ap-Af)!计?算?推?进?剂®的?燃?面?变?化ˉ规?律°并￠输?出?结®果?de=e1/numopen(10,file="export_burnS.dat")open(20,file="export_Ap.dat")do i=0,num,1if((i*de)<=r1) thenSa(i)=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(r1+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l-(r1-i*de)*Pi/n/l)*l*len Apa(i)=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-sin(epsilon*Pi/n)* cotan(thet/2)))+&2*n*(r+i*de)*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n)/l+&n*(r+i*de)**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2+n*(r1-i*de)**2*(thet/2+cotan(thet/2 )-Pi/2)/l**2)*l**2else if((i*de)>r1.and.(i*de)<=(l*sin(epsilon*Pi/n)/cos(thet/2)-r)) thenSa(i)=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(i*de+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l)*l*lenApa(i)=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-sin(epsilon*Pi/n)* cotan(thet/2)))+&2*n*(r+i*de)*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n)/l+&n*(r+i*de)**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2)*l**2else if((i*de)>(l*sin(epsilon*Pi/n)/cos(thet/2)-r).and.(i*de)<=e1) thenSa(i)=2*n*((1-epsilon)*Pi/n+(r+i*de)*(Pi/n+asin(l*sin(epsilon*Pi/n)/(i*de+r)))/l)*l *lenApa(i)=n*((1-epsilon)*Pi*(1+(r+i*de)/l)**2/n+sin(epsilon*Pi/n)*(sqrt((r+i*de)**2/l* *2-&sin(epsilon*Pi/n)**2)+cos(epsilon*Pi/n))+(r+i*de)**2*(epsilon*Pi/n+&asin(l*sin(epsilon*Pi/n)/(r+i*de)))/l**2)*l**2write(10,"(f10.5,2X,f15.5)") i*de,Sa(i)write(20,"(f10.5,2X,f15.5)") i*de,Apa(i)end doclose(10)close(20)!调獭?用?子哩?程²序´计?算?装痢?药?内¸弹獭?道台?曲¸线?callinternal_ballistics0(d,len,e1,n,thet,epsilon,r,r1,l,a,pn,Ctz,At,k,density_p,mpeff,Isp) !输?出?星?型³装痢?药?的?几?何?参?数簓open(30,file="export_star_geometry.dat")write(30,"(A6,f11.5,A2)") "mpeff=",mpeff,"Kg"write(30,"(A2,f8.5,A2)") "D=",d,"m"write(30,"(A2,f8.5,A2)") "L=",len,"m"write(30,"(A2,I2)") "n=",nwrite(30,"(A3,f8.5)") "θ¯=",thetwrite(30,"(A3,f7.5)") "ε?=",epsilonwrite(30,"(A2,f7.5,A2)") "r=",r,"m"write(30,"(A3,f7.5,A2)") "r1=",r1,"m"write(30,"(A3,f8.5,A2)") "e1=",e1,"m"write(30,"(A3,f10.5,A2)") "Ap=",Ap,"㎡O"write(30,"(A2,f10.5)") "J=",Jwrite(30,"(A3,f7.5)") "η?=",etawrite(30,"(A3,f10.5,A2)") "Af=",Af,"㎡O"write(30,"(A4,f10.5)") "η?f=",etafwrite(30,"(A2,f10.5)") "l=",lwrite(30,"(A3,f10.5)") "y0=",y0write(30,"(A3,f10.5)") "y1=",y1close(30)stopend!该?子哩?程²序´用?于 ?计?算?零?维?变?截?面?燃?烧?装痢?药?的?内¸弹獭?道台?!利?用?龙ⅷ?格?-库a塔t法ぁ?计?算?内¸弹獭?道台?曲¸线?!可°用?于 ?计?算?侵?蚀骸?燃?烧?效§应畖下?的?内¸弹獭?道台?曲¸线?!--------------------------------------------------------------------------!d外猘径?,len长¤度¯,e1平?均·肉╝厚?,n星?角?数簓,thet星?边?夹D角?,epsilon角?度¯系μ数簓!r过y度¯圆2弧?半?径?,r1星?角?圆2弧?半?径?和³l特?征¶尺?寸?,key表括?示?是?否?考?虑?侵?蚀骸?燃?烧?!density_p推?进?剂®的?密¹度¯,k比括?热¯?比括?ga系μ数簓,C特?征¶速·度¯,mpeff药?柱·质±量?!a速·度¯系μ数簓,pn压1强?指?数簓,At喉³部?面?积y,Pc燃?烧?室酣?的?设Θ?计?压1力 ,Isp表括?示?理え?论?比括?冲?!dt时骸?间?步?长¤,e燃?层?厚?度¯,time时骸?间?,Sa同?一?时骸?刻²轴®向´的?各¶节¸点?的?燃?面?,Apa通 ?气?面?积y!Poc各¶节¸点?的?压1力 ,Pe喷?管¹出?口¸压1力 ,P_I压1力 冲?量?,ep侵?蚀骸?比括?!P_av平?均·压1力 ,F_av平?均·推?力 ,I0总哩?冲?,Im重?量?比括?冲?Iv体?积y比括?冲?,Cf推?力 系μ数簓!--------------------------------------------------------------------------subroutineinternal_ballistics0(d,len,e1,n,thet,epsilon,r,r1,l,a,pn,Ctz,At,k,density_p,mpeff,Isp)implicit nonereal,parameter::Pi=3.14integer :: n,keyreal(kind=8),intent(in) :: d,len,e1,thet,epsilon,r,r1,lreal(kind=8),intent(in) :: density_p,k,Ctz,mpeff,Ispreal(kind=8) :: a,pn,At,dt,e,time,Sa,Apa,Poc,Pe,P_I,F real(kind=8) :: ep,P_av,F_av,I0,Im,Iv,Cf,gareal(kind=8) :: f1,f2,f3,f4Pe=101325.0P_I=0.0!药?柱·的?能¹量?特?性?参?量?读®入?open(3,file="inter_parameter.dat")read(3,*) keyread(3,*) dtread(3,*) Pocclose(3)!计?算?点?火e压1强?时骸?的?推?力 值μPoc=Poc*101325.0ga=sqrt(k)*(2/(k+1))**((k+1)/(k-1)/2)Cf=ga*sqrt(2*k*(1-(Pe/Poc)**((k-1)/k))/(k-1))F=Cf*Poc*Ate=0.0time=0.0open(10,file="export_Poc.dat")open(20,file="export_F.dat")write(10,"(f8.5,2X,f10.5)") time,Poc/101325.0write(20,"(f8.5,2X,f15.5)") time,Fdo while(.true.)!计?算?燃?面?和³通 ?气?面?积yif(e<=r1) thenSa=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(r1+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l-(r1-e)*Pi/n/l)*l*len Apa=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-&sin(epsilon*Pi/n)*cotan(thet/2)))+2*n*(r+e)*(sin(epsilon*Pi/n)/&sin(thet/2)+(1-epsilon)*Pi/n)/l+n*(r+e)**2*(Pi/n+Pi/2-thet/2-&cotan(thet/2))/l**2+n*(r1-e)**2*(thet/2+cotan(thet/2)-Pi/2)/l**2)*l**2else if(e>r1.and.e<=(l*sin(epsilon*Pi/n)/cos(thet/2)-r)) thenSa=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(e+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l)*l*lenApa=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-&sin(epsilon*Pi/n)*cotan(thet/2)))+2*n*(r+e)*(sin(epsilon*Pi/n)/&sin(thet/2)+(1-epsilon)*Pi/n)/l+n*(r+e)**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2)*l**2 else if(e>(l*sin(epsilon*Pi/n)/cos(thet/2)-r).and.e<=e1) thenSa=2*n*((1-epsilon)*Pi/n+(r+e)*(Pi/n+asin(l*sin(epsilon*Pi/n)/(e+r)))/l)*l*len Apa=n*((1-epsilon)*Pi*(1+(r+e)/l)**2/n+sin(epsilon*Pi/n)*(sqrt((r+e)**2/l**2-& sin(epsilon*Pi/n)**2)+cos(epsilon*Pi/n))+(r+e)**2*&(epsilon*Pi/n+asin(l*sin(epsilon*Pi/n)/(r+e)))/l**2)*l**2elseSa=0.0end if!计?算?侵?蚀骸?比括?if(key==1) thenif(Sa/Apa<=72.9) thenep=1.0elseep=1.3128-1.3249e-2*Sa/Apa+1.5527e-4*(Sa/Apa)**2-4.3868e-7*(Sa/Apa)**3end ifelseep=1.0end if!采°用?4阶¬的?龙ⅷ?格?库a塔t法ぁ?计?算?f1=ga**2*Ctz**2*(density_p*Sa*ep*a*Poc**pn-Poc*At/Ctz)/(Apa*len)f2=ga**2*Ctz**2*(density_p*Sa*ep*a*(Poc+dt*f1/2)**pn-(Poc+dt*f1/2)*At/Ctz)/(Apa*len )f3=ga**2*Ctz**2*(density_p*Sa*ep*a*(Poc+dt*f2/2)**pn-(Poc+dt*f2/2)*At/Ctz)/(Apa*len )f4=ga**2*Ctz**2*(density_p*Sa*ep*a*(Poc+dt*f3)**pn-(Poc+dt*f3)*At/Ctz)/(Apa*len) Poc=Poc+dt*(f1+2.0*f2+2.0*f3+f4)/6.0Cf=ga*sqrt(2*k*(1-(Pe/Poc)**((k-1)/k))/(k-1))F=Cf*Poc*Attime=time+dte=e+dt*ep*a*Poc**pnif(e>e1.and.Poc<101325.0) exitwrite(10,"(f8.5,2X,f10.5)") time,Poc/101325.0write(20,"(f8.5,2X,f15.5)") time,FP_I=P_I+Poc*dtend doclose(10)close(20)!循-环«结®束?P_av=P_I/timeCf=ga*sqrt(2*k*(1-(Pe/P_av)**((k-1)/k))/(k-1))F_av=Cf*P_av*AtI0=F_av*timeIm=I0/mpeff/9.81Iv=I0/(Apa*len)!输?出?固²体?火e箭y发ぁ?动ˉ机¸的?工¤作痢?参?数簓值μopen(30,file="export_propellant.dat")write(30,"('平?均·压1力 p_av=',f9.5)") p_av/101325.0write(30,"('平?均·推?力 F_av=',f10.5,1X,'KN')") F_av/1000.0 write(30,"('总哩?冲?I0=',f10.5,1X,'KN.s')") I0/1000.0write(30,"('重?量?比括?冲?Im=',f10.5,1X,'s')") Imwrite(30,"('体?积y比括?冲?Iv=',f10.5,1X,'KN.s/m3')") Iv/1000.0 write(30,"('工¤作痢?时骸?间?time=',f8.5,1X,'s')") timeclose(30)returnend subroutine。