一种考虑烧蚀-沉积的固体姿轨控发动机零维内弹道计算方法

1.本发明涉及姿轨控发动机性能方针领域，具体是一种考虑烧蚀-沉积的固体姿轨控发动机零维内弹道计算方法。


背景技术：

2.固体姿轨控发动机是现代动能武器的直接力控制系统，也是未来提升战术导弹武器系统作战性能的关键技术之一。固体姿轨控发动机通常由多喷管耦合工作改变推力矢量大小及方向，实现导弹的大推力快速变轨和姿态调节。发动机推力根据调节阀运动而实时改变，通过调整喉栓与喷管的相对位置直接影响发动机等效喉部面积，从而影响喷管质量流率和推力大小。发动机燃烧室压强是极其重要的设计参数，也是保证发动机稳定工作的必要条件，其变化规律直接决定了发动机的推力及工作时间。
3.目前，通常采用常规零维内弹道对燃烧室压强进行计算，不考虑等效喉部面积由于烧蚀/沉积出现的偏差。根据已知喉栓母线、喷管母线以及初始时刻喉栓位置，计算可得喉栓位移-等效喉部面积变化关系。通过对已知喉栓位移-等效喉部面积关系插值，得出各个阀门当前喉栓位移及所对应等效喉部面积，叠加不同时刻各喷管等效喉部面积可得等效喉部面积总和，代入零维内弹道基本方程中可得工作状态下燃烧室压强。
4.现有固体姿轨控发动机零维内弹道求解方法主要来源于固体发动机零维内弹道求解方法，计算发动机各种工作条件下燃烧室内压强随时间和空间的变化规律。在理想状态下，固体姿轨控发动机等效喉部面积可简单视为各阀门等效喉部面积之和，但实际工作过程中发动机喷管喉部气流复杂流动，传热和沉积、烧蚀的问题严重，采用常规零维内弹道对固体姿轨控发动机燃烧室压强进行求解，仅能得出理想状态下燃烧室压强情况，会造成燃烧室理论压强与实际压强偏差较大、推力控制不精确的问题。


技术实现要素：

5.针对上述现有技术中由于发动机燃温较高、燃料不纯净、内部材料不耐高温，在喉栓头部及喷管喉部出现烧蚀、沉积效应，造成发动机喷管型面改变，实际等效喉部面积与理想值出现偏差，影响发动机推力性能，造成固体姿轨控发动机推力调控不精确的问题，本发明提供一种考虑烧蚀-沉积的固体姿轨控发动机零维内弹道计算方法，能够进行已知烧蚀/沉积规律后的固体姿轨控发动机零维内弹道计算，获得更贴近工程实际情况的燃烧室压强，降低由于零维内弹道计算不精确造成的推力计算误差，易于控制机构对推力的精准调控，便于后期对发动机进行性能考核和数值模拟。
6.为实现上述目的，本发明提供一种考虑烧蚀-沉积的固体姿轨控发动机零维内弹道计算方法，包括如下步骤：
7.建立考虑烧蚀/沉积影响后的固体姿轨控发动机零维内弹道计算模型；
8.采用龙哥库塔法对固体姿轨控发动机零维内弹道计算模型进行迭代求解，得到考
虑烧蚀/沉积影响后对应燃烧室的压强；
9.其中，在每步龙格库塔法迭代前，基于喉栓位移及喷管、喉栓的烧蚀/沉积速率得到当前各阀门的等效喉部面积，并累加得到发动机的总等效喉部面积。
10.在另一个实施例中，所述考虑烧蚀/沉积影响后的固体姿轨控发动机零维内弹道计算模型为：
[0011][0012]
式中，pc是燃烧室压强，t是时间，γ是比热比k的函数，c
*
是特征速度，vc是燃烧室自由容积，ρ
p
是推进剂密度，ab是装药的燃面面积，a是燃速系数，n是压强指数，a
t
是各阀门等效喉部面积之和，e是药柱肉厚，r是燃气的气体常数，tf是定压燃烧温度，k是比热比，a
ti
是单阀等效喉部面积，xi是喉栓位移，r
t
、r
p
分别是喷管喉部和喉栓的烧蚀/沉积量。
[0013]
在另一个实施例中，采用嵌套二分法得到当前各阀门的等效喉部面积，在单次等效喉部面积求解过程中分为外循环层与内循环层，其中，外循环层在喉栓型面曲线上进行二分法搜索，内循环层在喷管型面曲线上进行二分法搜索。
[0014]
在另一个实施例中，所述外循环层的过程为：
[0015]
步骤1.1，获取当前喉栓位移位置，并根据烧蚀/沉积模型更新考虑烧蚀/沉积效应后的喉栓型面曲线，得到第一搜索区间；
[0016]
步骤1.2，选取第一搜索区间的左、右端点以及中点附近的任意两点作为第一控制点，并基于内循环层得到各第一控制点对应的等效喉部面积；
[0017]
步骤1.3，筛选出步骤1.2中等效喉部面积最小的第一控制点，并更新第一搜索区间，以该最小等效喉部面积对应的第一控制点左、右相邻的两个第一控制点作为第一搜索区间更新后的左、右端点；
[0018]
步骤1.4，判断第一搜索区间是否达到收敛条件，若是则输出步骤1.2中最小的等效喉部面积作为对应阀门的等效喉部面积，否则返回步骤1.2。
[0019]
在另一个实施例中，步骤1.3中，所述基于内循环层得到各第一控制点对应的等效喉部面积，具体为：
[0020]
对于每一第一控制点，都有：
[0021]
步骤2.1，根据烧蚀/沉积模型更新考虑烧蚀/沉积效应后的喷管型面曲线，得到第二搜索区间；
[0022]
步骤2.2，选取第二搜索区间的左、右端点以及中点附近的任意两点作为第二控制点，并获取第一控制点与各第二控制点之间的等效喉部面积；
[0023]
步骤2.3，筛选出步骤2.2中等效喉部面积最小的第二控制点，并更新第二搜索区间，以该最小等效喉部面积对应的第二控制点左、右相邻的两个第二控制点作为第二搜索区间更新后的左、右端点；
[0024]
步骤2.4，判断第二搜索区间是否达到收敛条件，若是则输出步骤2.2中最小的等效喉部面积作为对应第一控制点的等效喉部面积，否则返回步骤2.2。
[0025]
在另一个实施例中，步骤2.2中，所述获取第一控制点与各第二控制点之间的等效喉部面积，具体为：
[0026][0027]
式中，a
ti
为等效喉部面积，f1(m,h
t
)与f2(n,xi,h
p
)分别为考虑烧蚀/沉积后喷管母线方程与喉栓母线方程，m为喷管上对应第二控制点的横坐标，h
t
为喷管喉部表面烧蚀/沉积物厚度，n为喉栓上对应第一控制点的横坐标，xi为喉栓顶点当前位移，h
p
为喉栓表面烧蚀/沉积物厚度。
[0028]
与现有技术相比，本发明提供的一种考虑烧蚀-沉积的固体姿轨控发动机零维内弹道计算方法，具有如下有益技术效果：
[0029]
1、考虑发动机燃烧过程中由于烧蚀/沉积效应造成的等效喉部面积变化，从而实现固体姿轨控发动机工作过程燃烧室压强精确、快速计算；
[0030]
2、将喉栓与喷管两部分烧蚀/沉积速率分开考虑，可实现两部分沉积特性不同情况下的零维内弹道计算，为实现发动机推力精确计算奠定技术基础。
附图说明
[0031]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
[0032]
图1为本发明实施例中考虑烧蚀/沉积效应的零维内弹道计算方法流程图；
[0033]
图2为本发明实施例中嵌套二分法搜索等效喉部面积流程图；
[0034]
图3为本发明实施例中等效喉部面积计算坐标系示意图；
[0035]
图4为本发明实施例中示例的轨控发动机布局示意图；
[0036]
图5为本发明实施例中示例的喷管喉栓型面结构示意图；
[0037]
图6为本发明实施例中示例的各阀门喉栓位移变化曲线图；
[0038]
图7为本发明实施例中示例的各阀门理想等效喉部面积变化曲线图；
[0039]
图8为本发明实施例中示例的各阀门考虑沉积效应后等效喉部面积变化曲线图；
[0040]
图9为本发明实施例中示例的燃烧室压强曲线图。
[0041]
本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
[0042]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0043]
需要说明，本发明实施例中所有方向性指示(诸如上、下、左、右、前、后
……
)仅用于解释在某一特定姿态(如附图所示)下各部件之间的相对位置关系、运动情况等，如果该特定姿态发生改变时，则该方向性指示也相应地随之改变。
[0044]
另外，在本发明中如涉及“第一”、“第二”等的描述仅用于描述目的，而不能理解为指示或暗示其相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
[0045]
在本发明中，除非另有明确的规定和限定，术语“连接”、“固定”等应做广义理解，例如，“固定”可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接，还可以是物理连接或无线通信连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系，除非另有明确的限定。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
[0046]
另外，本发明各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
[0047]
在固体姿轨控发动机工作过程中，由于发动机燃温较高、燃料不纯净、内部材料不耐高温等因素，在喉栓头部及喷管喉部会出现烧蚀、沉积效应，造成发动机喷管型面改变，实际等效喉部面积与理想值出现偏差，影响发动机推力性能，造成固体姿轨控发动机推力调控不精确的问题。基于此，本实施例以固体发动机内弹道学为基础，公开了一种高效、快速的考虑烧蚀-沉积的固体姿轨控发动机零维内弹道计算方法，该方法在考虑由烧蚀、沉积对等效喉部面积的影响后，采用嵌套二分法快速、准确地按时间步长更新实际等效喉部面积，实现了已知构型、已知固体姿轨控发动机阀门烧蚀/沉积规律、已知喉栓运动方案后的发动机燃烧室压强求解。按该方法进行已知烧蚀/沉积规律后的固体姿轨控发动机零维内弹道计算，可获得更贴近工程实际情况的燃烧室压强，降低由于零维内弹道计算不精确造成的推力计算误差，易于控制机构对推力的精准调控，便于后期对发动机进行性能考核和数值模拟，解决了发动机零维内弹道精确快速计算问题。相较于常规零维内弹道方法，采用本方法进行燃烧室压强计算，可获得更贴近实际工作过程且精度更高的结果，为固体姿轨控发动机推力精准调控奠定技术基础。
[0048]
本实施例首先将已知构型的固体姿轨控发动机喉栓、喷管两部分烧蚀/沉积效率进行定义，得到烧蚀/沉积模型。其次在已知喉栓运动方案后，将考虑烧蚀/沉积效应后的母线参数化表示。然后选用嵌套二分法对各时刻不同阀门的等效喉部面积进行快速精确求解，得出当前时刻考虑烧蚀/沉积后实际等效喉部面积总和，并代入零维内弹道计算模型中进行迭代计算。最终实现从初始时刻到最终时刻按时间步长迭代的烧蚀/沉积后零维内弹道计算。
[0049]
本实施例中，虑烧蚀-沉积的固体姿轨控发动机零维内弹道计算方法的具体实施过程如下：
[0050]
首先获取给定的姿轨控发动机构型参数、材料物性参数与计算设置参数。其中，姿轨控发动机构型参数包括装药燃面和质量特性随燃烧厚度的变化规律、喷管和喉栓型面参数；材料物性参数包括推进剂参数、喷管烧蚀/沉积模型和喉栓烧蚀/沉积模型；计算设置参数包括内弹道计算时间步长、喉栓运动方案。
[0051]
再建立考虑烧蚀/沉积影响后的固体姿轨控发动机零维内弹道计算模型，为：
[0052][0053]
式中，pc是燃烧室压强，t是时间，γ是比热比k的函数，c
*
是特征速度，vc是燃烧室自由容积，ρ
p
是推进剂密度，ab是装药的燃面面积，a是燃速系数，n是压强指数，a
t
是各阀门等效喉部面积之和，e是药柱肉厚，r是燃气的气体常数，tf是定压燃烧温度，k是比热比，a
ti
是单阀等效喉部面积，xi是喉栓位移，r
t
、r
p
分别是喷管喉部和喉栓的烧蚀/沉积量。等效喉部计算时应综合喉栓位移及喷管、喉栓的烧蚀/沉积速率进行单步的等效喉部面积计算。出现烧蚀情况时，实际等效喉部面积大于理论等效喉部面积；出现沉积情况时，实际等效喉部面积小于理论等效喉部面积。
[0054]
根据上述考虑烧蚀/沉积影响后的固体姿轨控发动机零维内弹道计算模型，参考图1，固体姿轨控发动机零维内弹道计算过程为：
[0055]
第1步：输入发动机及喷管参数、推进剂性能参数、燃面参数及喉栓运动方案；
[0056]
第2步：给定计算所需初值及常量；
[0057]
第3步：确定燃烧时间划分(给出时间步长)；
[0058]
第4步：初始燃烧室压强计算；
[0059]
第5步：确定当前时刻的喉栓、喷管烧蚀/沉积量，并计算等效喉部面积；
[0060]
第6步：燃速计算；
[0061]
第7步：燃去肉厚计算；
[0062]
第8步：对已有燃面数据进行插值得出对应燃面大小；
[0063]
第9步：当前喉栓运动方案考虑烧蚀-沉积影响后对应燃烧室压强计算；
[0064]
第10步：判断是否燃烧到药柱末端，若是则结束，否则返回第5步。
[0065]
本实施例中，采用四阶龙哥库塔法对固体姿轨控发动机零维内弹道计算模型进行迭代求解，得到考虑烧蚀/沉积影响后对应燃烧室的压强。其中，燃烧室初始压强由计算设置参数得到，装药燃面面积ab(e)与燃烧室自由容积vc(e)直接由燃面计算程序计算得到不同燃烧厚度e下的燃面面积和自由容积插值得到，不同压强和等效喉部面积条件下的阀门和喷管烧蚀/沉积速率h(pc,a
ti
)和g(pc,a
ti
)根据实验数据实测获得，各阀等效喉部面积a
ti
＝f(xi,r
t
,r
p
)通过步骤(4)所示的嵌套二分法获得，在每步龙格库塔法迭代前，基于喉栓位移及喷管、喉栓的烧蚀/沉积速率得到当前各阀门的等效喉部面积，并累加得到发动机的总等效喉部面积。
[0066]
采用嵌套二分法得到任意时刻不同阀门喉栓位移后、烧蚀/沉积后，快速计算其所对应的等效喉部面积总和。由于在喉栓在喷管运动过程中等效喉部面积存在唯一极值点，选用二分法对喉栓与喉部间最小通道面积进行搜索，能够大幅提高求解效率、减少计算消耗。等效喉部面积求解问题中由于计算的是喉栓和喷管母线之间的最小通道面积，应分别考虑两者的求解区间，因此本实施例采用双层嵌套的二分法进行循环搜索求解，在单次等效喉部面积求解时分为外循环层与内循环层，其中，外循环层在喉栓型面曲线上进行二分法搜索，内循环层在喷管型面曲线上进行二分法搜索。
[0067]
参考图2，外循环层的过程为：
[0068]
步骤1.1，获取当前喉栓位移位置，并根据烧蚀/沉积模型更新考虑烧蚀/沉积效应后的喉栓型面曲线，得到第一搜索区间；
[0069]
步骤1.2，选取第一搜索区间的左、右端点以及中点附近的任意两点作为第一控制点，并基于内循环层得到各第一控制点对应的等效喉部面积；
[0070]
步骤1.3，筛选出步骤1.2中等效喉部面积最小的第一控制点，并更新第一搜索区间，以该最小等效喉部面积对应的第一控制点左、右相邻的两个第一控制点作为第一搜索区间更新后的左、右端点；
[0071]
步骤1.4，判断第一搜索区间是否达到收敛条件，若是则输出步骤1.2中最小的等效喉部面积作为对应阀门的等效喉部面积，否则返回步骤1.2。其中，判断第一搜索区间是否达到收敛条件的方式为达到最大外循环迭代次数或第一搜索区间当前的长度小于预设的判断精度。
[0072]
在步骤1.3中，基于内循环层得到各第一控制点对应的等效喉部面积，对于每一第一控制点，都有：
[0073]
步骤2.1，读入外循环层中当前喉栓型面曲线的第一控制点，根据烧蚀/沉积模型更新考虑烧蚀/沉积效应后的喷管型面曲线，得到第二搜索区间；
[0074]
步骤2.2，选取第二搜索区间的左、右端点以及中点附近的任意两点作为第二控制点，并获取第一控制点与各第二控制点之间的等效喉部面积，其实施过程为：
[0075]
已知喉栓和喷管的型面参数，假设考虑其烧蚀/沉积规律后喉栓和喷管的型面函数，为：
[0076][0077]
式中，y1＝f1(m,h
t
)和y2＝f2(n,xi,h
p
)分别为考虑烧蚀/沉积后喷管母线方程与喉栓母线方程；m为喷管上任意一点横坐标，h
t
为喷管喉部表面烧蚀/沉积物厚度，m
l
、mr为喷管母线左右端点，即第一搜索区间的初始左、右端点；n为喉栓上任意一点横坐标，h
p
为喉栓表面烧蚀/沉积物厚度，n
l
、xi为喉栓母线左右端点，即第二搜索区间的初始左、右端点，xi也为喉栓顶点当前位移。
[0078]
计算喉栓的有效行程区间和喉栓不同运动位置对应的等效喉部面积的坐标系如图3所示。喷管喉部与喉栓表面沉积物厚度与沉积速率的关系为：
[0079]ht
＝r
t
δt
[0080]hp
＝r
p
δt
[0081]
式中，r
t
为喷管喉部表面烧蚀/沉积物沉积速率，r
p
为喉栓表面烧蚀/沉积物沉积速率，δt为时间步长。
[0082]
针对单阀而言，随着喉栓从左向右进动，等效喉部面积随喉栓的移动逐渐减小，喉栓与喉部完全闭合时达到最小值。当某喷管喉栓运动至位移xi处，计算其喉部面积a
ti
通常在喷管上任取一点(m,y1)、在喉栓上任取一点(n,y2)后计算等效喉部面积，遍历所有点，过程中求解出的最小面积为该阀等效喉部面积。遍历阀门位移xi的有效运动区间后，过程中所有位移值及所对应求解出的等效喉部面积，即为单阀喉栓位移-等效喉部面积变化关系。本实施例中，采用通过二分法来进行等效喉部面积搜索，能够大幅提高求解效率、减少计算消耗。等效喉部面积计算过程为：
[0083][0084]
式中，a
ti
为等效喉部面积，上式中的m即为喷管上对应第二控制点的横坐标，n即为喉栓上对应第一控制点的横坐标；
[0085]
步骤2.3，筛选出步骤2.2中等效喉部面积最小的第二控制点，并更新第二搜索区间，以该最小等效喉部面积对应的第二控制点左、右相邻的两个第二控制点作为第二搜索区间更新后的左、右端点；
[0086]
步骤2.4，判断第二搜索区间是否达到收敛条件，若是则输出步骤2.2中最小的等效喉部面积作为对应第一控制点的等效喉部面积，否则返回步骤2.2。其中，判断第二搜索区间是否达到收敛条件的方式为达到最大内循环迭代次数或第二搜索区间当前的长度小
于预设的判断精度。
[0087]
下面以某四喷管轨控发动机时序实验设计为例，给出一个具体案例进行方法验证，对本实施例中考虑烧蚀-沉积的固体姿轨控发动机零维内弹道计算方法作出进一步的说明。
[0088]
该实验设计中轨控发动机布局如图4所示，轨控发动机可产生沿本体oy轴与oz轴的轨控力，阀门按顺时针方向编号为阀门1、2、3、4。喷管喉栓型面如图5所示，其构型参数见表1，发动机工作参数如表2设置。
[0089]
表1喷管喉栓构型参数
[0090][0091][0092]
表2发动机工作参数表
[0093][0094]
该计算方法的目的是考虑固体姿轨控发动机工作过程中的烧蚀/沉积效应后，实现燃烧室压强的快速、精确计算，降低燃烧室压强理论计算值与实际值误差。针对该案例，应用考虑烧蚀/沉积效应的固体姿轨控发动机零维内弹道计算方法计算工作过程中燃烧室压强，具体步骤如下：
[0095]
步骤1、给出如图6所示喉栓运动方案，其对应理想等效喉部面积变化如图7所示。设定烧蚀/沉积模型如式所示，将阀门喉栓与喉部沉积规律分开考虑，以式喉栓表面沉积厚度计算公式与式喷管表面沉积厚度计算公式更新当前母线方程。在定义的工作时间内按照时间步长进行循环求解，已知0s时刻作为初始时刻，按照时间步长0.01s进行叠加，计算至最终时刻10s。
[0096]
hd＝rδt＝aa
t
(t)mδt
[0097]
其中，hd是烧蚀/沉积物的厚度，r是烧蚀/沉积速率，δt是时间步长，a
t
是理想等效喉部面积，a是烧蚀/沉积项系数，m是等效喉部面积变化影响烧蚀/沉积项变化指数。
[0098][0099][0100]
其中，h
pi_d
是各喉栓表面附着沉积物的厚度，h
ti_d
是各喷管表面附着沉积物的厚度，r
pi_d
是喉栓沉积速率，r
ti_d
是喷管沉积速率，δt是时间步长，a
ti
是各阀门理想等效喉部面积，pc是燃烧室理想压强，a1,a2分别是喉栓沉积项系数和喷管沉积项系数，m1,m2分别是等效喉部面积变化影响喉栓沉积项变化指数和等效喉部面积变化影响喷管沉积项变化指数。
[0101]
步骤2、采用嵌套二分法对考虑烧蚀/沉积效应后的喷管进行等效喉部面积计算，得出考虑沉积效应的等效喉部面积变化曲线，如图8所示。对比图7、图8可以看出，由于考虑沉积物堵塞对喉部面积有效变化范围的消耗，考虑沉积效应后的等效喉部面积变化区间逐渐减小，而理想情况下的等效喉部面积变化区间始终不变，考虑烧蚀/沉积效应的等效喉部面积变化更符合实际工程经验。
[0102]
步骤3、根据计算得出的考虑烧蚀/沉积效应后的等效喉部面积计算零维内弹道计算，得到考虑沉积后的燃烧室压强，理想压强是不考虑烧蚀/沉积效应对等效喉部面积的影响计算得出的燃烧室压强，如图9所示。由图9可以明显看出理想压强与考虑沉积效应后的燃烧室压强有较大差异，若采用常规零维内弹道计算方法与发动机实际情况差异较大，采用本实施例方法计算得出的燃烧室压强更能支撑后期对发动机性能研究工作的开展。
[0103]
以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。