航天飞行动力学课程设计——飞船再入质点弹道日期：2022-04-27航天飞行动力学课程设计 0——飞船再入质点弹道 01.题目重述 (1)1)假设：12)标称轨迹制导 12.背景分析 (2)3.数值求解方法 (2)1)地球以及大气模型22)再入初始数据 23)线性插值方法 24)积分方法-四阶龙格库塔 25)蒙特卡洛打靶随机数生成24.分析过程 (3)1)求解ODE获取基准弹道 32)给定偏差量求解ODE获取制导弹道弹道35.结果分析 (3)1)基准弹道情况 32)100次打靶结果分析56.C++程序结构及主要代码 (6)1)头文件62)Cpp文件63)函数声明 74)函数定义 81. 题目重述1) 假设：● 考虑地球旋转影响。
● 地球看成质量均匀分布的圆球，质心在球心。
● 把飞行器看成质点，应用瞬时平衡假设。
2222sin cos sin cos cos cos sin cos (sin cos cos sin cos )1cos ()cos 2cos sin cos (cos cos sin cos sin )1sin cos sin tan 2cos e e e drV dt d V dt r d V dt r dV D g r dt d V L g V r dt V r d L V dt V r γθγψφφγψγωφγφγφψγσγωφψωφγφγψφψσγψφγ====--+-⎡⎤=+-+++⎢⎦⎣⎡=+-⎢⎣2(1)(tan cos cos sin )sin sin cos cos e e r V ωωγψφφψφφγ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎤⎪-+⎥⎪⎦⎩ 上述动力学方程组中，有6个状态变量：[,,,,,]r V θφγψ。
各状态变量的意义为：r ：地球球心到飞行器质心的距离；λ：经度；φ：纬度；V ：相对地球速度；γ：速度倾角；ψ：速度方位角，0ψ=表示正北方向，从正北顺时针旋转为正。
e ω为地球旋转角速度；,D L 分别为阻力加速度和升力加速度，可由下式给出：2211(,)(,)(2)22ref D ref L D V S C Ma L V S C Ma mmραρα==,D LC C 分别为飞行器的阻力系数和升力系数，它们是攻角α和马赫数的函数；ref S为飞行器参考面积；ρ为大气密度。
首先按照配平攻角飞行，得到基准弹道。
2) 标称轨迹制导倾侧角指令(/)cos /cL D L Dσ=0(/)(/)(/)c L D L D L D =+∆，其中0(/)L D 为基准弹道升阻比，取为0.28；(/)L D ∆为与以速度为自变量的基准弹道偏差引起的升阻比，由下式计算：1234(/)x L D k n k R k h k R ∆=∆+∆+∆+∆x n ∆为切向过载偏差，R ∆为航程偏差。
1234,,,k k k k 为系数，通过试验法自行确定。
倾侧角指令在轴向过载大于0.5的时候开始输出，在轴向过载小于0.5时，采用开环制导的方式,即常数10度。
2. 背景分析制导背景：该飞船再入使用弹道-升力式再入，通过配置质心的方法，使航天器进入大气层时产生一定升力，其质心不配置在再入航天器的中心轴线上，而偏离中心轴线一小段距离，同时质心在压心之前，故航天器使用配平攻角飞行，此升力一般不大于阻力的一半，即升阻比小于0.5，其精度比弹道式更优良，外形为简单的旋成体，在一定范围内可以控制航天器的着陆点位置，其最大过载也大大小于弹道式再入时的最大过载。
动力学背景：以配平攻角飞行时，空气动力R 通过航天器的压心和质心，且再入航天器为旋成体，其压心在再入航天器的几何纵轴上，侧滑角为0，攻角小于0.3. 数值求解方法1) 地球以及大气模型重力模型g =g 0(1−(h Re)2)其中g 0=9. 9.80665m/s 2,Re 为地球半径6378.137Km 。
地球自转角速度w ie =7.292e-05rad/s;大气密度模型(10km~120km)选用《航天飞行动力学》P.294~P.295d 的USSA1976标准大气表的拟合值。
声速公式按如下公式计算a =20.0468√T式中的温度（K ）使用大气密度模型进行插值。
2) 再入初始数据=9500m=238.ref S0120km h =θ=0-4; φ=-015 =07600Vγ=03-ψ=0533) 线性插值方法利用Ma 对平衡攻角，阻力系数，升力系数，升阻比进行一阶线性插值。
此外，调取基准弹道偏差量时，也需要进行以速度V 为自变量的线性插值。
4) 积分方法-四阶龙格库塔11234612122322243(22)(,)(,)(,)(,)h i i i i h hi i h hi i i i x x K K K K K f t x K f t x K K f t x K K f t h x hK +=++++==++=++=++5) 蒙特卡洛打靶随机数生成多元线性同余算法生成随机数：本算法选用的线性组合系数a=[269,113,17]，全为质数时性能更加；选用m=16384，默认x0=[91,5,13]。
4. 分析过程1) 求解ODE 获取基准弹道根据6维ODE 方程组：进行4阶龙格库塔积分可求出6个状态变量：[,,,,,]r V θφγψ。
于是可以绘制相关变量关于时间的变化曲线和基准弹道曲线。
方程组中的D 和L 是与Ma 有关的变量，可通过插值得到，(/)cos /cL D L Dσ=，0(/)(/)(/)c L D L D L D =+∆(/)L D ∆按0计算，0(/)L D 取为0.28。
2) 给定偏差量求解ODE 获取制导弹道弹道对飞船加上制导和统计偏差后，1234(/)x L D k n k R k h k R ∆=∆+∆+∆+∆，通过试验调节k1,k2,k3,k4，来控制制导精度，使落点偏差尽可能小，同时避免飞船发散，统计偏差用随机数生成，最后确定落点偏差的均值和置信区间。
本次打靶，对于K 值的初步的设计结果如下k1 = -2e-1, k2 = -5e-7, k3 = -1e-3; k4 = -2e-5;由于本次大作业时间仓促，仅仅只对于这四个值做了初步的设计，只要保证弹道不发散而已。
然而实际使用过程中，必须判断K 值取的好坏。
标准是，能否把任意范围内的拉偏量调整回基准弹道附近，使最终的脱靶量最小，落点精度最高。
5. 结果分析1) 基准弹道情况绘制基准弹道结果曲线，高度-时间, 速度-时间,迎角-时间, 动压-时间,过载-时间，弹道倾角-时间，航程-时间，阻力加速度-时间，倾侧角曲线；基准弹道为S形曲线，在一段时间内飞船平飞，高度变化不大，航程变化率逐渐降低。
速度随高度降低而下降，动压和过载的变化规律相似，在较大的范围内变化。
以上变量均是波动变化的，可以分析出切向过载主要是和阻力加速度相关的。
分析：攻角随高度的降低缓慢增加，而在降低到某一高度时，攻角快速降低。
弹道倾角前期变化不大，而在降低到某一高度时随快速降低。
倾侧角在切向过载小于0.5时输出常值10°，降低到某一高度时变化范围极大，与末制导段复杂的气动特性有关。
2)256次打靶结果分析对终端落点偏差进行蒙特卡洛打靶分析，确定落点偏差的均值和置信区间。
绘制蒙特卡洛打靶的相关曲线。
在加入制导和统计偏差后，进行了256次打靶实验，每一次打靶的弹道如下：对于这256次打靶落点，将其投影到LLH坐标系，可以直观地看到打靶结果的落点散布：可以发现落点基本在以标准弹道落点附近的30km之内。
落点偏差的均方差估计值(25934,22097)m 落点偏差在置信区间（-5930.12，+5930.12）内的概率为95%分析：在k1,k2,k3,k4取值适当的情况下，能得到教理想的弹道，在不同的偏差下，精度基本符合要求。
6.C++程序结构及主要代码本程序创建了三个个类C_RK4和C_linPol以及C_MultiLinMod分别用来实现4阶Runge-Kutta积分、线性插值以及多元线性随机数生成。
使用的时候包含这三个类的头文件和实现文件（cpp文件），就可以创建这三个类的实例；在main函数调用对象的成员函数就可以实现相应的计算以及输出。
这样的对象化编程可以避免直接调用函数时出现太多的形参表，进而简洁高效，有逻辑地实现设定被积函数、初始值、步长、触发条件，以及使用不同方式进行积分和输出结果。
此外，本算法利用C++的重载特性，可以实现不同边界条件下的积分。
调用对象Reentry的不同成员函数，可以实现不同的数值积分算法。
以下介绍程序的结构：1)头文件●C_RK4.h 类声明●C_linPol.h 类声明●C_MultiLinMod.h 类声明●Pch.h 函数以及全局变量声明2)Cpp文件●C_RK4.cpp 4阶Runge-Kutta积分的定义文件●C_linPol.cpp 线性插值定义文件●C_MultiLinMod.cpp 多元线性同余法随机数生成器●Pch.cpp 主要定义全局变量、定义被积函数和它所调用的函数●Main.cpp 执行文件由于所使用的三个类都是通用的文件，以下只对于本次工程项目直接要使用的文件做一说明。
3)#ifndef PCH_H#define PCH_H#include<iostream>//#include <iomanip>#include <fstream>#include<sstream>#include<vector>#include <cstdlib>#include <cmath>#include "C_RK4.h"#include "C_linPol.h"#include "CMulitLinMod.h"// 静态常量—全局变量声明const int neq = 7;static const double g0 = 9.80665, w_ie = 7.2921e-05, Re = 6378137, d2r = 180.0 / 3.141592654;// 函数—求解ODE所需要的函数float atmosphere_density(float h, float *Temperature);float gravity(float h);// 重力模型float soundSpeed(float T);// 声速模型float balanceAOA(double L_D_qS[3], float V, float h);// 配平攻角求解升力阻力void bias2StdTrj(double *y, double biasLpD[]);// 求解相对于基准弹道的偏差double gudiance(float LpD, float N_x);// 制导方程（无偏差）double gudiance(float LpD, float N_x, double t, double *y); 制导方程（打靶使用）void reentry(double t, double *y, double *yDerivative);// ODE方程组float above10km(double *y, float dt);// 判别终止条件int readStdTrjFile(double *stdV, double *stdGamma, double *stdnx, double *stdR);// 读取基准弹道void prt(int num, double *R);// 测试读取值是否正确void MonteCarlo(C_RK4 &trj, const int times);// 打靶并保存结果#endif //PCH_H4)函数定义// pch.cpp: 与预编译标头对应的源文件；编译成功所必需的#include "pch.h"// 未扰动的初始值const float m = 9500, S_ref = 23.8,V0 = 7600, gamma0 = -3.0;static float dm = 0, dL = 0, dD = 0, dTheta = 0, dV = 0;// 在这里定义的作用区域更小，只限于pch.cpp。