航天飞行动力学课程设计-飞船再入质点弹道数值计算Word文档下载推荐.doc
《航天飞行动力学课程设计-飞船再入质点弹道数值计算Word文档下载推荐.doc》由会员分享，可在线阅读，更多相关《航天飞行动力学课程设计-飞船再入质点弹道数值计算Word文档下载推荐.doc（16页珍藏版）》请在冰点文库上搜索。
资料
1.题目重述
1）假设：
l考虑地球旋转影响。
l地球看成质量均匀分布的圆球，质心在球心。
l把飞行器看成质点，应用瞬时平衡假设。
上述动力学方程组中，有6个状态变量：
。
各状态变量的意义为：
：
地球球心到飞行器质心的距离；
经度；
纬度；
相对地球速度；
速度倾角；
速度方位角，表示正北方向，从正北顺时针旋转为正。
为地球旋转角速度；
分别为阻力加速度和升力加速度，可由下式给出：
分别为飞行器的阻力系数和升力系数，它们是攻角和马赫数的函数；
为飞行器参考面积；
为大气密度。
首先按照配平攻角飞行，得到基准弹道。
2）标称轨迹制导
倾侧角指令
，
其中为基准弹道升阻比，取为0.28；
为与以速度为自变量的基准弹道偏差引起的升阻比，由下式计算：
为切向过载偏差，为航程偏差。
为系数，通过试验法自行确定。
倾侧角指令在轴向过载大于0.5的时候开始输出，在轴向过载小于0.5时，采用开环制导的方式,即常数10度。
2.背景分析
制导背景：
该飞船再入使用弹道-升力式再入，通过配置质心的方法，使航天器进入大气层时产生一定升力，其质心不配置在再入航天器的中心轴线上，而偏离中心轴线一小段距离，同时质心在压心之前，故航天器使用配平攻角飞行，此升力一般不大于阻力的一半，即升阻比小于0.5，其精度比弹道式更优良，外形为简单的旋成体，在一定范围内可以控制航天器的着陆点位置，其最大过载也大大小于弹道式再入时的最大过载。
动力学背景：
以配平攻角飞行时，空气动力R通过航天器的压心和质心，且再入航天器为旋成体，其压心在再入航天器的几何纵轴上，侧滑角为0，攻角小于0.
3.数值求解方法
1）地球以及大气模型
重力模型
其中g0=9.9.80665m/s2,Re为地球半径6378.137Km。
地球自转角速度wie=7.292e-05rad/s;
大气密度模型（10km~120km）选用《航天飞行动力学》P.294~P.295d的USSA1976标准大气表的拟合值。
声速公式按如下公式计算
式中的温度（K）使用大气密度模型进行插值。
2）再入初始数据
;
3）线性插值方法
利用Ma对平衡攻角，阻力系数，升力系数，升阻比进行一阶线性插值。
此外，调取基准弹道偏差量时，也需要进行以速度V为自变量的线性插值。
4）积分方法-四阶龙格库塔
5）蒙特卡洛打靶随机数生成
多元线性同余算法生成随机数：
本算法选用的线性组合系数a=[269,113,17]，全为质数时性能更加；
选用m=16384，默认x0=[91,5,13]。
4.分析过程
1）求解ODE获取基准弹道
根据6维ODE方程组：
进行4阶龙格库塔积分可求出6个状态变量：
于是可以绘制相关变量关于时间的变化曲线和基准弹道曲线。
方程组中的D和L是与Ma有关的变量，可通过插值得到，
按0计算，取为0.28。
2）给定偏差量求解ODE获取制导弹道弹道
对飞船加上制导和统计偏差后，，通过试验调节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之内。
落点偏差的均方差估计值
落点偏差在置信区间（-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）头文件
lC_RK4.h 类声明
lC_linPol.h 类声明
lC_MultiLinMod.h 类声明
lPch.h 函数以及全局变量声明
2）Cpp文件
lC_RK4.cpp 4阶Runge-Kutta积分的定义文件
lC_linPol.cpp 线性插值定义文件
lC_MultiLinMod.cpp 多元线性同余法随机数生成器
lPch.cpp 主要定义全局变量、定义被积函数和它所调用的函数
lMain.cpp 执行文件
由于所使用的三个类都是通用的文件，以下只对于本次工程项目直接要使用的文件做一说明。
变量名的定义如下：
3）函数声明
#ifndefPCH_H
#definePCH_H
#include<
iostream>
//#include<
iomanip>
#include<
fstream>
sstream>
vector>
cstdlib>
cmath>
#include"
C_RK4.h"
C_linPol.h"
CMulitLinMod.h"
//静态常量—全局变量声明
constintneq=7;
staticconstdoubleg0=9.80665,w_ie=7.2921e-05,Re=6378137,d2r=180.0/3.141592654;
//函数—求解ODE所需要的函数
floatatmosphere_density（floath,float*Temperature）;
floatgravity（floath）;
// 重力模型
floatsoundSpeed（floatT）;
// 声速模型
floatbalanceAOA（doubleL_D_qS[3],floatV,floath）;
// 配平攻角求解升力阻力
voidbias2StdTrj（double*y,doublebiasLpD[]）;
// 求解相对于基准弹道的偏差
doublegudiance（floatLpD,floatN_x）;
// 制导方程（无偏差）
doublegudiance（floatLpD,floatN_x,doublet,double*y）;
制导方程（打靶使用）
voidreentry（doublet,double*y,double*yDerivative）;
// ODE方程组
floatabove10km（double*y,floatdt）;
// 判别终止条件
intreadStdTrjFile（double*stdV,double*stdGamma,double*stdnx,double*stdR）;
// 读取基准弹道
voidprt（intnum,double*R）;
// 测试读取值是否正确
voidMonteCarlo（C_RK4&
trj,constinttimes）;
// 打靶并保存结果
#endif//PCH_H
4）函数定义
//pch.cpp:
与预编译标头对应的源文件；
编译成功所必需的
pch.h"
//未扰动的初始值
constfloatm=9500,S_ref=23.8,V0=7600,gamma0=-3.0;
staticfloatdm=0,dL=0,dD=0,dTheta=0,dV=0;
//在这里定义的作用区域更小，只限于pch.cpp。
而在pch.h里边定义，相当于直接定义了全局变量，作用区域包括本文件和main函数
inttnuit=10;
//配平攻角一维插值源数据
staticdoubleMa_index[15]={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.2414,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_linPolinterp_CD（Ma_index,Cd_mat）;
C_linPolinterp_CL（Ma_index,Cl_mat）;
C_linPolinterp_AOA（Ma_index,AOA_mat）;
C_linPolinterp_std_gamma;
C_linPolinterp_std_nx;
C_linPolinterp_std_R;
//所调用的函数。
floatatmosphere_density（floath,float*Temperature）{
/*测试函数如下
floatP,T;
for（inti=0;
i<
120;
i++）{
P=atmosphere_density（i,&
T）;
std:
:
cout<
<
'
'
<
P<
T<
\n'
}*/
doubleH,W,T,P,R0=6371.393;
doubleps=1.225;
H=h/（1+h/R0）;
if（h>
=0&
h<
=11.0191）{
W=1-H/44.3308;
T=288.15*W;
P=pow（W,4.2559）*ps;
}
elseif（h>
11.0191&
h<
=20.0631）
{
W=exp（（14.9647-H）/6.3416）;
T=216.650;
P=0.15898*W*ps;
20.0631&
=32.1619）
W=1+（H-24.9021）/221.552;
T=221.552*W;
P=3.2722*0.01*pow（W,-35.1629）*ps;
32.1619&
=47.3501）
W=1+（H-39.7499）/89.4107;
T=250.350*W;
P=3.2618*0.001*pow（W,-13.2011）*ps;
47.3501&
=51.4125）
W=exp（（48.6252-H）/7.9223）;
T=270.650;
P=9.4920*pow（10,-4）*W*ps;
51.4125&
=71.8020）
W=1-（H-59.4390）/88.2218;
T=247.021*W;
P=2.5280*pow（10,-4）*pow（W,11.2011）*ps;
71.8020&
=86.0000）
W=1-（H-78.0303）/100.2950;
T=200.59*W;
P=1.7632*pow（10,-5）*pow（W,16.0816）*ps;
86.0000&
=91.0000）
W=exp（（87.2848-H）/5.47）;
T=186.87;
P=3.6411*pow（10,-6）*W*ps;
else
P=ps*exp（-h/7.320）/3;
}*Temperature=T;
returnP;
}
floatgravity（floath）{
returng0*（1-（h/Re）*（h/Re））;
floatsoundSpeed（floatT）{
return20.0468*sqrt（T）;
floatbalanceAOA（doubleD_L_qS[3],floatV,floath）{
floatrho,T,qS,Mach;
rho=atmosphere_density（h/1000.0,&
qS=0.5*V*V*rho*S_ref;
Mach=V/soundSpeed（T）;
doubleCl,Cd;
Cl=interp_CL.intern（Mach）,Cd=interp_CD.intern（Mach）;
D_L_qS[0]=Cd*（1+dL）*qS/（m+dm）;
D_L_qS[1]=Cl*（1+dD）*qS/（m+dm）;
D_L_qS[2]=qS;
returninterp_AOA.intern（Mach）;
//攻角的返回值单位为°
voidbias2StdTrj（double*y,doublebiasLpD[]）{
doublenx=y[7],nx_Std,gamma_Std,gamma=y[4];
floatV=y[3],R=y[6],R_Std;
nx_Std=interp_std_nx.intern（（double）V）;
gamma_Std= interp_std_gamma.intern（（double）V）;
R_Std=interp_std_R.intern（（double）V）;
biasLpD[0]=nx-nx_Std;
biasLpD[1]=R-R_Std;
biasLpD[2]=（cos（gamma）-cos（gamma_Std））*V;
biasLpD[3]=（sin（gamma）-sin（gamma_Std））*V;
doublegudiance（floatLpD,floatN_x）{
doubleLpD_c=0.28;
if（N_x>
0.5）returnLpD_c/LpD;
elsereturn0.984807753012;
//=cos（10°
）
doublegudiance（floatLpD,floatN_x,doublet,double*y）{
doubleLpD_c=0.28,biasLpD[4];
if（t>
118）{
intstop;
if（abs（N_x）>
0.5）{
floatk1,k2,k3,k4;
k1=-5.136e-1,k3=-1.503e-3,k2=-5.218e-6;
k4=-2.179e-5;
bias2StdTrj（y,biasLpD）;
LpD_c+=（k1*biasLpD[0]+k2*biasLpD[1]+k3*biasLpD[2]+k4*biasLpD[3]）;
returnLpD_c/LpD;
voidreentry（doublet,double*y,double*yDerivative）
{
/*常微分方程
变量按次序为 运动学：
r,lambda,phi,
动力学：
V,gamma,psi
*/
doubleh=y[0],lambda=y[1],phi=y[2],V=y[3],gamma=y[4],psi=y[5],alpha,D_L_qS[3]={0,0,0};
floatr=h+Re,g=gravity（h）;
alpha=balanceAOA（D_L_qS,V,h）;
doubleD=D_L_qS[0],L=D_L_qS[1];
floatqS=D_L_qS[2],N_x=D*cos（alpha/d2r）-L*sin（alpha/d2r）;
N_x=N_x/g;
//轴向力
//纵向制导
doublecs=gudiance（L/D,N_x,t,y）,ss=sqrt（1-cs*cs）;
doublesin_gamma=sin（gamma）,cos_gamma=cos（gamma）,cos_phi=cos（phi）,sin_phi=sin（phi）,cos_psi=cos（psi）,sin_psi=sin（psi）;
tnuit）{
tnuit+=10;
//微分方程组
yDerivative[0]=V*sin_gamma;
//高度
yDerivative[6]=V*cos_gamma;
//水平射程
yDerivative[1]=yDerivative[6]*sin_psi/（r*cos_phi）;
yDerivative[2]=yDerivative[6]*cos_psi/r;
yDerivative[3]=-D-g*sin_gamma+w_ie*w_ie*r*cos_phi*（sin_gamma*cos_phi-cos_gamma*sin_phi*cos_psi）;
yDerivative[4]=（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;
yDerivative[5]=（L*ss/cos_gamma+V*V/r*cos_gamma*sin_psi*tan（phi）-2*w_ie*V*（tan（gamma）*cos_psi*cos_p
升级会员 