航天飞行动力学课程设计飞船再入质点弹道数值计算Word格式文档下载.docx

1、倾侧角指令，其中 为基准弹道升阻比，取为0.28；为与以速度为自变量的基准弹道偏差引起的升阻比，由下式计算：为切向过载偏差，为航程偏差。为系数，通过试验法自行确定。倾侧角指令在轴向过载大于0.5的时候开始输出，在轴向过载小于0.5时，采用开环制导的方式,即常数10度。背景分析制导背景：该飞船再入使用弹道-升力式再入，通过配置质心的方法，使航天器进入大气层时产生一定升力，其质心不配置在再入航天器的中心轴线上，而偏离中心轴线一小段距离，同时质心在压心之前，故航天器使用配平攻角飞行，此升力一般不大于阻力的一半，即升阻比小于0.5，其精度比弹道式更优良，外形为简单的旋成体，在一定范围内可以控制航天器的

2、着陆点位置，其最大过载也大大小于弹道式再入时的最大过载。动力学背景：以配平攻角飞行时，空气动力R通过航天器的压心和质心，且再入航天器为旋成体，其压心在再入航天器的几何纵轴上，侧滑角为0，攻角小于0.数值求解方法3） 地球以及大气模型重力模型其中g0=9. 9.80665m/s2,Re为地球半径6378.137Km。地球自转角速度wie=7.292e-05rad/s;大气密度模型（10km120km）选用航天飞行动力学P.294P.295d的USSA1976标准大气表的拟合值。声速公式按如下公式计算式中的温度（K）使用大气密度模型进行插值。4） 再入初始数据 ;5） 线性插值方法利用Ma对平衡攻

3、角，阻力系数，升力系数，升阻比进行一阶线性插值。此外，调取基准弹道偏差量时，也需要进行以速度V为自变量的线性插值。6） 积分方法-四阶龙格库塔7） 蒙特卡洛打靶随机数生成多元线性同余算法生成随机数：本算法选用的线性组合系数a=269,113,17，全为质数时性能更加；选用m=16384，默认x0=91,5,13。分析过程8） 求解ODE获取基准弹道根据6维ODE方程组：进行4阶龙格库塔积分可求出6个状态变量：于是可以绘制相关变量关于时间的变化曲线和基准弹道曲线。方程组中的D和L是与Ma有关的变量，可通过插值得到，按0计算，取为0.28。9） 给定偏差量求解ODE获取制导弹道弹道对飞船加上制导和

4、统计偏差后，通过试验调节k1,k2,k3,k4，来控制制导精度，使落点偏差尽可能小，同时避免飞船发散，统计偏差用随机数生成，最后确定落点偏差的均值和置信区间。本次打靶，对于K值的初步的设计结果如下k1 = -2e-1, k2 = -5e-7, k3 = -1e-3; k4 = -2e-5;由于本次大作业时间仓促，仅仅只对于这四个值做了初步的设计，只要保证弹道不发散而已。然而实际使用过程中，必须判断K值取的好坏。标准是，能否把任意范围内的拉偏量调整回基准弹道附近，使最终的脱靶量最小，落点精度最高。结果分析10） 基准弹道情况绘制基准弹道结果曲线，高度-时间, 速度-时间,迎角-时间, 动压-时间

5、,过载-时间，弹道倾角-时间，航程-时间，阻力加速度-时间，倾侧角曲线；基准弹道为S形曲线，在一段时间内飞船平飞，高度变化不大，航程变化率逐渐降低。速度随高度降低而下降，动压和过载的变化规律相似，在较大的范围内变化。以上变量均是波动变化的，可以分析出切向过载主要是和阻力加速度相关的。分析：攻角随高度的降低缓慢增加，而在降低到某一高度时，攻角快速降低。弹道倾角前期变化不大，而在降低到某一高度时随快速降低。倾侧角在切向过载小于0.5时输出常值10，降低到某一高度时变化范围极大，与末制导段复杂的气动特性有关。11） 256次打靶结果分析对终端落点偏差进行蒙特卡洛打靶分析，确定落点偏差的均值和置信区间

6、。绘制蒙特卡洛打靶的相关曲线。在加入制导和统计偏差后，进行了256次打靶实验，每一次打靶的弹道如下：对于这256次打靶落点，将其投影到LLH坐标系，可以直观地看到打靶结果的落点散布：可以发现落点基本在以标准弹道落点附近的30km之内。落点偏差的均方差估计值落点偏差在置信区间（-5930.12，+5930.12）内的概率为95%在k1,k2,k3,k4取值适当的情况下，能得到教理想的弹道，在不同的偏差下，精度基本符合要求。C+程序结构及主要代码本程序创建了三个个类C_RK4和C_linPol以及C_MultiLinMod分别用来实现4阶Runge-Kutta积分、线性插值以及多元线性随机数生成。

7、使用的时候包含这三个类的头文件和实现文件（cpp文件），就可以创建这三个类的实例；在main函数调用对象的成员函数就可以实现相应的计算以及输出。这样的对象化编程可以避免直接调用函数时出现太多的形参表，进而简洁高效，有逻辑地实现设定被积函数、初始值、步长、触发条件，以及使用不同方式进行积分和输出结果。此外，本算法利用C+的重载特性，可以实现不同边界条件下的积分。调用对象Reentry的不同成员函数，可以实现不同的数值积分算法。以下介绍程序的结构：12） 头文件 C_RK4.h 类声明 C_linPol.h 类声明 C_MultiLinMod.h 类声明 Pch.h 函数以及全局变量声明13） C

8、pp文件 C_RK4.cpp 4阶Runge-Kutta积分的定义文件 C_linPol.cpp 线性插值定义文件 C_MultiLinMod.cpp 多元线性同余法随机数生成器 Pch.cpp 主要定义全局变量、定义被积函数和它所调用的函数 Main.cpp 执行文件由于所使用的三个类都是通用的文件，以下只对于本次工程项目直接要使用的文件做一说明。变量名的定义如下：14） 函数声明#ifndef PCH_H#define PCH_H#include/#include #include sstreamvectorcstdlibcmath#include C_RK4.hC_linPol.hCMu

9、litLinMod.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_qS3, f

10、loat 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, floa

11、t 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_H15） 函数定义/ pch.cpp: 与预编译标头对应的源文件；编译成功所必需的pch.h/ 未扰动的初始值const float m = 9500, S_ref = 23.8,V0

12、 = 7600, gamma0 = -3.0;static float dm = 0, dL = 0, dD = 0, dTheta = 0, dV = 0;/ 在这里定义的作用区域更小，只限于pch.cpp。而在pch.h里边定义，相当于直接定义了全局变量，作用区域包括本文件和main函数int tnuit = 10;/ 配平攻角一维插值源数据static double Ma_index15 = 0.5,0.7,0.9,1.1,1.2,1.5,2,2.4,3,4,6,10,18,25,32.2 ,Cd_mat = 0.86879,0.92765,1.0116,1.1366,1.717,1.2

13、414,1.2282,1.1808,1.1525,1.1444,1.1651,1.1886,1.2307,1.2446,1.2507 ,Cl_mat = 0.32106,0.26019,0.35029,0.51974,0.54683,0.57913,0.53645,0.5186,0.51709,0.50069,0.47066,0.46326,0.44825,0.4376,0.43513 ,AOA_mat = 17.43,17.54,22.36,27.97,29.22,29.67,29.45,29.2,28.65,27.2,26.68,25.94,24.18,23.68,23.22 ;C_lin

14、Pol interp_CD（Ma_index, Cd_mat）;C_linPol interp_CL（Ma_index, Cl_mat）;C_linPol interp_AOA（Ma_index, AOA_mat）;C_linPol interp_std_gamma;C_linPol interp_std_nx;C_linPol interp_std_R;/ 所调用的函数。float atmosphere_density（float h,float *Temperature） /* 测试函数如下 float P,T; for （int i = 0; i 120; i+） P=atmospher

15、e_density（i, &T）; std:cout , P T = 0 & h 11.0191&h 0.5）return LpD_c / LpD; else return 0.984807753012;/ =cos（10）double gudiance（float LpD, float N_x, double t, double *y） double LpD_c = 0.28,biasLpD4; if （t 118） int stop; if （abs（N_x） 0.5） float k1, k2, k3, k4; k1 = -5.136e-1, k3 = -1.503e-3, k2 = -

16、5.218e-6; k4 = -2.179e-5; bias2StdTrj（y, biasLpD）; LpD_c += （k1*biasLpD0 + k2 * biasLpD1 + k3 * biasLpD2 + k4 * biasLpD3）; return LpD_c / LpD;void reentry（double t, double *y, double *yDerivative） /*常微分方程 变量按次序为 运动学：r,lambda,phi, 动力学：V,gamma,psi */ double h = y0, lambda = y1, phi = y2, V = y3, gamma

17、 = y4, psi = y5,alpha, D_L_qS3 = 0,0,0 ; float r = h + Re, g = gravity（h）; alpha=balanceAOA（D_L_qS, V, h）; double D = D_L_qS0, L = D_L_qS1; float qS = D_L_qS2,N_x=D*cos（alpha/d2r）-L*sin（alpha/d2r）; N_x = N_x / g;/ 轴向力 / 纵向制导 double cs=gudiance（L/D,N_x,t,y）, ss=sqrt（1-cs*cs）; double sin_gamma = sin（g

18、amma）, cos_gamma = cos（gamma）, cos_phi = cos（phi）, sin_phi = sin（phi）, cos_psi = cos（psi）, sin_psi = sin（psi）; tnuit） tnuit += 10; / 微分方程组 yDerivative0 = V * sin_gamma;/ 高度 yDerivative6 = V * cos_gamma;/ 水平射程 yDerivative1 = yDerivative6 *sin_psi / （r*cos_phi）; yDerivative2 = yDerivative6 *cos_psi /

19、r; yDerivative3 = -D - g * sin_gamma + w_ie * w_ie*r*cos_phi*（sin_gamma*cos_phi - cos_gamma*sin_phi*cos_psi）; yDerivative4 = （L*cs + （V*V / r - g）*cos_gamma + 2 * w_ie*V*cos_phi*sin_psi + w_ie * w_ie*r*cos_phi*（cos_gamma*cos_phi + sin_gamma*cos_psi*sin_phi） / V; yDerivative5 = （L*ss / cos_gamma + V

20、* V / r * cos_gamma*sin_psi*tan（phi） - 2 * w_ie*V*（tan（gamma）*cos_psi*cos_phi - sin_phi） + w_ie * w_ie*r*sin_psi*sin_phi*cos_phi / cos_gamma） / V; /yDerivative7 = 1;/ 时间 ，便与调试 /yDerivative8 = yDerivative3; / 这一部分用于保存计算过程中的中间变量 yDerivative7 = yDerivative3/g;/ 切向过载 yDerivative8 = qS;/ 动压 yDerivative9

21、= D;/ 阻力 yDerivative10 = acos（cs）;/ sigma yDerivative11 = alpha;/ 攻角float above10km（double *y,float dt） / 终止条件 / 函数指针仍要调用这个函数 double Vy=y3 * sin（y4）; /if （y1 0.253355350897）return 0.001; /if （abs（y4） 1.5&abs（y5） 1.5）return 0.001; if （Vy*dt+ y0） 10000.0） return （10000.0 -y0） / Vy; return 10000;int readStdTrjFile（double *stdV, double *stdGamma, double *stdnx, doub