本发明涉及可视化仿真技术领域，具体而言，涉及一种可视化弹道模型的设计方法。
背景技术：
现有的飞行器弹道设计一般采用基于C、C++、Matlab等编程语言的自编程序，技术人员在进行弹道设计时面临的是成千上万行程序代码，不容易上手，效率较低，人机界面不友好、可读性较差，同时不同设计人员编制的程序各不相同，通用性较差。
技术实现要素：
本发明实施方式的目的在于提供一种可视化弹道模型的设计方法，其操作简单，上手容易，弹道模型的设计效率高，通用性好。
本发明的实施方式是这样实现的：
本发明的实施方式提供了一种可视化弹道模型的设计方法，包括以下步骤：
步骤一：定义模型单元
利用面向对象的语言将弹道模型中的每个方程都单独定义成一个模型单元；
步骤二：定义关系类
利用面向对象的语言定义不同模型单元之间的关系类，每个关系类定义了两个模型单元之间的输入输出关系，每个关系类定义的两个模型单元中的一个模型单元为输入模型单元，另一个模型单元为输出模型单元；
步骤三：排列模型单元和关系类的运算次序；
步骤四：给定积分步长和积分算法，根据步骤三中确定的运算次序进行运算，运算完成输出计算结果。
进一步地，步骤三中所述的排列模型单元和关系类的运算次序包含以下过程：
S1：新建空队列；
S2：将步骤一中定义的不需要输入变量的模型单元存入到空队列中，如果有多个模型单元都不需要输入，则按照建立模型单元的先后顺序将其存入空队列中；
S3：根据步骤二中定义的关系类判断剩余模型单元的所有输入模型单元是否已在空队列中，如果不在空队列中则将其存入空队列中，直至没有模型单元可以存入为止，如果所有的模型单元都已经存入空队列，则空队列中模型单元的先后次序就是积分运算次序。
进一步地，在S3中，当没有模型单元可以存入空队列中时，如果有剩余的输入模型单元没有存入空队列中，则说明所构建的弹道模型存在代数环，不能求解，给出错误信息并返回步骤一重新定义该模型单元。
进一步地，步骤四中所述的根据步骤三中确定的运算次序进行运算包含以下过程：
S4：给定时间变量t、积分步长step、积分初始时间t0和积分结束时间tf，并令t＝t0；
S5：根据运算次序运算所有的模型单元；
S6：令t＝t+step，当t＜tf时，返回S5；当t≥tf时，运算结束并输出计算结果。
进一步地，积分步长为0.01-0.02s，积分算法为欧拉法或龙格库塔积分算法。
优选地，积分步长为0.01s。
进一步地，所述面向对象的语言包括C++、C#或Java。
进一步地，步骤一中所述的弹道模型中的方程包括质心运动学方程、质心动力学方程、姿态运动学方程和姿态动力学方程。
本发明的实施方式具有以下有益效果：
本发明实施方式提供的可视化弹道模型的设计方法，其操作简单，上手容易，设计人员通过简单的选取和输入即可完成弹道模型搭建，并且可以将公知的和一些重复使用的模型单元统一封装，并预留通用的输入输出接口，进而实现“搭积木”式弹道建模，大大提高了弹道模型的设计效率；根据所定义的各模型单元之间的关系类，能够自动排列各模型单元和关系类的运算次序，给定积分步长和积分算法，从而达到正确求解弹道模型的目的。
附图说明
为了更清楚地说明本发明实施方式的技术方案，下面将对实施方式中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施方式，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1为本发明实施方式提供的可视化弹道模型的设计方法的流程图；
图2为关系类的定义示意图；
图3为运算次序排列示意图；
图4为发动机推力模型单元(T)的计算示意图；
图5为质量模型单元(A)的计算示意图；
图6为攻角、侧滑角模型单元(B)的计算示意图；
图7为大气模型单元(C)的计算示意图；
图8为动压模型单元(D)的计算示意图；
图9为马赫数模型单元(E)的计算示意图；
图10为质心、转动惯量特性模型单元(F)的计算示意图；
图11为气动模型单元(G)的计算示意图；
图12为发射系下的加速度模型单元(H)的计算示意图；
图13为绕质心转动模型单元(I)的计算示意图；
图14为姿态运动学模型单元(J)的计算示意图；
图15为关系类的定义示意图；
图16为可视化弹道设计模型及弹道计算结果示意图；
图17为飞行器高度曲线示意图。
图标：101-第一模型单元；102-第二模型单元；103-第三模型单元；104-模型单元a；105-模型单元b；106-模型单元c；107-模型单元d；108-模型单元e；109-模型单元f。
具体实施方式
为了加深对本发明的理解，下面结合实施方式对本发明作进一步说明，该实施方式仅用于解释本发明，并不构成对本发明保护范围的限定。
参考图1所示，本发明实施方式提供的一种可视化弹道模型的设计方法，包括以下步骤：
步骤一：定义模型单元
利用面向对象的语言将弹道模型中的每个方程都单独定义成一个模型单元。
需要说明的是，本实施方式中，面向对象的语言包括C++、C#或Java等语言。
步骤二：定义关系类
利用面向对象的语言定义不同模型单元之间的关系类，每个关系类定义了两个模型单元之间的输入输出关系，每个关系类定义的两个模型单元中的一个模型单元为输入模型单元，另一个模型单元为输出模型单元。
参考图2所示，第一模型单元101的输入为x1和x2，第二模型单元102的输出为x1，第三模型单元103的输出为x2，则可以建立两个关系类，这两个关系类分别为第二模型单元102到第一模型单元101的关系类和第三模型单元103到第一模型单元101的关系类，其中，第一模型单元101为输入模型单元，第二模型单元102和第三模型单元103均为输出模型单元。
步骤三：排列模型单元和关系类的运算次序。具体包含以下过程：
S1：新建空队列，用于存放模型单元；
S2：将步骤一中定义的不需要输入变量的模型单元存入到空队列中，如果有多个模型单元都不需要输入，则按照建立模型单元的先后顺序将其存入空队列中；
S3：根据步骤二中定义的关系类判断剩余模型单元的所有输入模型单元是否已在空队列中，如果不在空队列中则将其存入空队列中，直至没有模型单元可以存入为止，如果所有的模型单元都已经存入空队列，则空队列中模型单元存入的先后次序就是积分运算次序；如果有剩余的输入模型单元没有存入空队列中，则说明所构建的弹道模型存在代数环，不能求解，给出错误信息并返回步骤一。
参考图3所示，依次存入空队列的模型单元是第一次存入模型单元a104，第二次存入模型单元b105，第三次放入模型单元c106，第四次放入模型单元d107，第五次放入模型单元e108和模型单元f109,则积分运算次序为先运算模型单元a104，然后依次运算模型单元b105、模型单元c106、模型单元d107、模型单元e108和模型单元f109。
步骤四：给定积分步长和积分算法，根据步骤三中确定的运算次序进行运算，运算完成输出计算结果。具体包含以下过程：
S4：给定时间变量t、积分步长step、积分初始时间t0和积分结束时间tf，并令t＝t0；
S5：根据运算次序运算所有的模型单元；
S6：令t＝t+step，当t＜tf时，返回S5；当t≥tf时，运算结束并输出计算结果。
需要说明的是，本实施方式中，积分步长为0.01-0.02s，优选积分步长为0.01s；积分算法为欧拉法或龙格库塔积分算法。
还需要说明的是，本实施方式中，弹道模型方程包括质心运动学方程、质心动力学方程、姿态运动学方程、姿态动力学方程、速度方程、速度方向方程、角度方程、角速度方程和高度方程等等。
下面通过算例进行说明：
已知某飞行器起飞前质量为M0，飞行过程中推进剂秒耗量为Is，发动机推力为Fp，飞行初始位置在发射坐标系中的值为X0＝(x0,y0,z0)，初始速度在发射坐标系中的值为V0＝(Vx0,Vy0,Vz0)，假设地球为不旋转的圆球模型。则弹道建模和计算过程如下：
一、定义模型单元
(1)发动机推力模型单元(T)，其求解方程为：
式中，F(t)和G(t)分别为发动机推力和秒耗量随时间的函数，一般为拟合曲线或插值表。
参考图4所示，向发动机推力模型单元(T)输入时间t，即可求解得到发动机推力Fp和推进剂秒耗量Is。
(2)质量模型单元(A)，其求解方程为：
式中，M为飞行器的质量。
参考图5所示，向质量模型单元(A)输入推进剂秒耗量Is，即可求解得到飞行器的质量M。
(3)攻角、侧滑角模型单元(B)，其求解方程为：
首先，根据当前的速度求取速度倾角θ和航迹偏航角σ：
式中，vx为飞行器速度在发射系x轴上的投影，vy为飞行器速度在发射系y轴上的投影，vz为飞行器速度在发射系z轴上的投影。
根据当前的三个姿态欧拉角(俯仰角)、ψ(偏航角)、γ(滚转角)和速度倾角θ、航迹偏航角σ求取攻角α和侧滑角β：
参考图6所示，向攻角、侧滑角模型单元(B)输入ψ、γ、Vx、Vy和Vz，即可求解得到攻角α和侧滑角β。
(4)大气模型单元(C)，其求解方程为：
根据飞行高度求取大气密度ρ和音速sonic：
式中，y为飞行器的飞行高度，f(y)和g(y)分别为大气密度随飞行高度变化的函数表达式和音速随飞行高度变化的函数，一般为拟合曲线或插值表。
参考图7所示，向大气模型单元(C)输入飞行器的高度y，即可求解得到大气密度ρ和音速sonic。
(5)动压模型单元(D)，其求解方程为：
式中，Q为飞行器飞行动压。
参考图8所示，向动压模型单元(D)输入ρ、Vx、Vy和Vz，即可求解得到Q。
(6)马赫数模型单元(E)，其求解方程为：
式中，Ma为飞行器飞行马赫数。
参考图9所示，向马赫数模型单元(E)输入音速sonic、Vx、Vy和Vz，即可求解得到Ma。
(7)质心、转动惯量特性模型单元(F)，其求解方程为：
以质量为插值变量，利用质量特性插值表插值得到质心和转动惯量：
式中，Xc、Yc、Zc分别为飞行器质心在弹体坐标系下三个方向的位置；Jx、Jy、Jz分别为飞行器绕弹体坐标系三个坐标轴的转动惯量；h1(M)、h2(M)、h3(M)、h4(M)、h5(M)和h6(M)分别为飞行器质心位置三分量和转动惯量三分量随质量变化的函数，一般为拟合曲线或插值表。
参考图10所示，向质心、转动惯量特性模型单元(F)输入飞行器的质量M，即可求解得到Xc、Yc、Zc、Jx、Jy和Jz。
(8)气动模型单元(G)，其求解方程为：
首先根据马赫数、攻角、侧滑角等输入插值得到气动力系数和气动力矩系数CA、CN、CZ、CMX、CMY、CMZ：
式中，i1(Ma，α，β)、i2(Ma，α，β)……i6(Ma，α，β)分别为气动六分量系数随马赫数、攻角、侧滑角变化的函数，一般为拟合曲线或插值表。
然后根据动压、质心位置及气动力和气动力矩系数求取气动力和气动力矩：
参考图11所示，向气动模型单元(G)输入α、β、Ma、Q、Xc，即可求解得到Rx1、Ry1、Rz1、Mx1、My1和Mz1。
(9)发射系下的加速度模型单元(H)，其求解方程为：
首先计算弹体系下三个方向的视加速度：
将弹体系下的视加速度变换到发射系下：
求取发射系下的加速度：
式中，Tvx、Tvy、Tvz分别为发射系下三个方向的加速度。
参考图12所示，向发射系下的加速度模型单元(H)输入Fp、Rx1、Ry1、Rz1、M、ψ和γ，即，可求解得到Tvx、Tvy和Tvz。
(10)绕质心转动模型单元(I)，其求解方程为：
参考图13所示，向绕质心转动模型单元(I)输入Mx1、My1、Mz1、Jx、Jy、Jz中，即可求解得到ωx、ωy和ωz。
(11)姿态运动学模型单元(J)，其求解方程为：
参考图14所示，向姿态运动学模型单元(J)输入ωx、ωy和ωz，即可求解得到发射系到弹体系的三个欧拉角ψ和γ。
二、定义关系类
梳理第一步中各模型单元的输入输出变量及各模型单元之间的相对关系，并将其定义为关系类，关系类的定义示意图入图15所示。本实施方式中定义了：发动机推力模型单元(T)到质量模型单元(A)的关系类，发动机推力模型单元(T)到发射系下的加速度模型单元(H)的关系类，攻角、侧滑角模型单元(B)到气动模型单元(G)的关系类，大气模型单元(C)到动压模型单元(D)的关系类，大气模型单元(C)到马赫数模型单元(E)的关系类，质量模型单元(A)到质心、转动惯量特性模型单元(F)的关系类，质量模型单元(A)到加速度模型单元(H)的关系类，动压模型单元(D)到气动模型单元(G)的关系类，马赫数模型单元(E)到气动模型单元(G)的关系类，质心、转动惯量特性模型单元(F)到绕质心转动模型单元(I)的关系类，质心、转动惯量特性模型单元(F)到气动模型单元(G)的关系类，气动模型单元(G)到加速度模型单元(H)的关系类，气动模型单元(G)到绕质心转动模型单元(I)的关系类。
三、排列模型单元和关系类的运算次序
根据第二步定义的关系类，可以得到运算次序如下：首先运算发动机推力模型单元(T)、攻角、侧滑角模型单元(B)、大气模型单元(C)、姿态运动学模型单元(J)，其次运算质量模型单元(A)、动压模型单元(D)、马赫数模型单元(E)，然后运算质心、转动惯量特性模型单元(F)，再次是运算气动模型单元(G)，最后是运算加速度模型单元(H)和绕质心转动模型单元(I)。
四、给定积分步长(例如0.01s)和积分算法(例如四阶龙格库塔积分算法)，根据第三步中确定的运算次序进行积分运算，运算完成输出计算结果，即可得到弹道计算结果。
五、可视化弹道设计模型及弹道计算结果如图16所示。图17为飞行器飞行高度的结果曲线。
以上所述仅为本发明的优选实施方式而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。