航天飞行动力学课程设计飞船再入质点弹道数值计算.docx
- 文档编号:10776237
- 上传时间:2023-02-22
- 格式:DOCX
- 页数:25
- 大小:344.67KB
航天飞行动力学课程设计飞船再入质点弹道数值计算.docx
《航天飞行动力学课程设计飞船再入质点弹道数值计算.docx》由会员分享,可在线阅读,更多相关《航天飞行动力学课程设计飞船再入质点弹道数值计算.docx(25页珍藏版)》请在冰豆网上搜索。
航天飞行动力学课程设计飞船再入质点弹道数值计算
航天飞行动力学课程设计
——飞船再入质点弹道
班级
02011502
姓名
XXX
XXXX
学号
**********
**********
日期:
2020-02-01
题目重述
1)假设:
●考虑地球旋转影响。
●地球看成质量均匀分布的圆球,质心在球心。
●把飞行器看成质点,应用瞬时平衡假设。
上述动力学方程组中,有6个状态变量:
。
各状态变量的意义为:
:
地球球心到飞行器质心的距离;
:
经度;
:
纬度;
:
相对地球速度;
:
速度倾角;
:
速度方位角,
表示正北方向,从正北顺时针旋转为正。
为地球旋转角速度;
分别为阻力加速度和升力加速度,可由下式给出:
分别为飞行器的阻力系数和升力系数,它们是攻角
和马赫数的函数;
为飞行器参考面积;
为大气密度。
首先按照配平攻角飞行,得到基准弹道。
2)标称轨迹制导
倾侧角指令
,
其中
为基准弹道升阻比,取为0.28;
为与以速度为自变量的基准弹道偏差引起的升阻比,由下式计算:
为切向过载偏差,
为航程偏差。
为系数,通过试验法自行确定。
倾侧角指令在轴向过载大于0.5的时候开始输出,在轴向过载小于0.5时,采用开环制导的方式,即常数10度。
背景分析
制导背景:
该飞船再入使用弹道-升力式再入,通过配置质心的方法,使航天器进入大气层时产生一定升力,其质心不配置在再入航天器的中心轴线上,而偏离中心轴线一小段距离,同时质心在压心之前,故航天器使用配平攻角飞行,此升力一般不大于阻力的一半,即升阻比小于0.5,其精度比弹道式更优良,外形为简单的旋成体,在一定范围内可以控制航天器的着陆点位置,其最大过载也大大小于弹道式再入时的最大过载。
动力学背景:
以配平攻角飞行时,空气动力R通过航天器的压心和质心,且再入航天器为旋成体,其压心在再入航天器的几何纵轴上,侧滑角为0,攻角小于0.
数值求解方法
3)地球以及大气模型
重力模型
其中g0=9.9.80665m/s2,Re为地球半径6378.137Km。
地球自转角速度wie=7.292e-05rad/s;
大气密度模型(10km~120km)选用《航天飞行动力学》P.294~P.295d的USSA1976标准大气表的拟合值。
声速公式按如下公式计算
式中的温度(K)使用大气密度模型进行插值。
4)再入初始数据
;
5)线性插值方法
利用Ma对平衡攻角,阻力系数,升力系数,升阻比进行一阶线性插值。
此外,调取基准弹道偏差量时,也需要进行以速度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获取制导弹道弹道
对飞船加上制导和统计偏差后,
,通过试验调节k1,k2,k3,k4,来控制制导精度,使落点偏差尽可能小,同时避免飞船发散,统计偏差用随机数生成,最后确定落点偏差的均值和置信区间。
本次打靶,对于K值的初步的设计结果如下
k1=-2e-1,k2=-5e-7,k3=-1e-3;k4=-2e-5;
由于本次大作业时间仓促,仅仅只对于这四个值做了初步的设计,只要保证弹道不发散而已。
然而实际使用过程中,必须判断K值取的好坏。
标准是,能否把任意范围内的拉偏量调整回基准弹道附近,使最终的脱靶量最小,落点精度最高。
结果分析
10)基准弹道情况
绘制基准弹道结果曲线,高度-时间,速度-时间,迎角-时间,动压-时间,过载-时间,弹道倾角-时间,航程-时间,阻力加速度-时间,倾侧角曲线;
基准弹道为S形曲线,在一段时间内飞船平飞,高度变化不大,航程变化率逐渐降低。
速度随高度降低而下降,动压和过载的变化规律相似,在较大的范围内变化。
以上变量均是波动变化的,可以分析出切向过载主要是和阻力加速度相关的。
分析:
攻角随高度的降低缓慢增加,而在降低到某一高度时,攻角快速降低。
弹道倾角前期变化不大,而在降低到某一高度时随快速降低。
倾侧角在切向过载小于0.5时输出常值10°,降低到某一高度时变化范围极大,与末制导段复杂的气动特性有关。
11)256次打靶结果分析
对终端落点偏差进行蒙特卡洛打靶分析,确定落点偏差的均值和置信区间。
绘制蒙特卡洛打靶的相关曲线。
在加入制导和统计偏差后,进行了256次打靶实验,每一次打靶的弹道如下:
对于这256次打靶落点,将其投影到LLH坐标系,可以直观地看到打靶结果的落点散布:
可以发现落点基本在以标准弹道落点附近的30km之内。
落点偏差的均方差估计值
落点偏差在置信区间(-5930.12,+5930.12)内的概率为95%
分析:
在k1,k2,k3,k4取值适当的情况下,能得到教理想的弹道,在不同的偏差下,精度基本符合要求。
C++程序结构及主要代码
本程序创建了三个个类C_RK4和C_linPol以及C_MultiLinMod分别用来实现4阶Runge-Kutta积分、线性插值以及多元线性随机数生成。
使用的时候包含这三个类的头文件和实现文件(cpp文件),就可以创建这三个类的实例;在main函数调用对象的成员函数就可以实现相应的计算以及输出。
这样的对象化编程可以避免直接调用函数时出现太多的形参表,进而简洁高效,有逻辑地实现设定被积函数、初始值、步长、触发条件,以及使用不同方式进行积分和输出结果。
此外,本算法利用C++的重载特性,可以实现不同边界条件下的积分。
调用对象Reentry的不同成员函数,可以实现不同的数值积分算法。
以下介绍程序的结构:
12)头文件
●C_RK4.h类声明
●C_linPol.h类声明
●C_MultiLinMod.h类声明
●Pch.h函数以及全局变量声明
13)Cpp文件
●C_RK4.cpp4阶Runge-Kutta积分的定义文件
●C_linPol.cpp线性插值定义文件
●C_MultiLinMod.cpp多元线性同余法随机数生成器
●Pch.cpp主要定义全局变量、定义被积函数和它所调用的函数
●Main.cpp执行文件
由于所使用的三个类都是通用的文件,以下只对于本次工程项目直接要使用的文件做一说明。
变量名的定义如下:
14)函数声明
#ifndefPCH_H
#definePCH_H
#include
//#include
#include
#include
#include
#include
#include
#include"C_RK4.h"
#include"C_linPol.h"
#include"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
15)函数定义
//pch.cpp:
与预编译标头对应的源文件;编译成功所必需的
#include"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<
}*/
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;
}
elseif(h>20.0631&h<=32.1619)
{
W=1+(H-24.9021)/221.552;
T=221.552*W;
P=3.2722*0.01*pow(W,-35.1629)*ps;
}
elseif(h>32.1619&h<=47.3501)
{
W=1+(H-39.7499)/89.4107;
T=250.350*W;
P=3.2618*0.001*pow(W,-13.2011)*ps;
}
elseif(h>47.3501&h<=51.4125)
{
W=exp((48.6252-H)/7.9223);
T=270.650;
P=9.4920*pow(10,-4)*W*ps;
}
elseif(h>51.4125&h<=71.8020)
{
W=1-(H-59.4390)/88.2218;
T=247.021*W;
P=2.5280*pow(10,-4)*pow(W,11.2011)*ps;
}
elseif(h>71.8020&h<=86.0000)
{
W=1-(H-78.0303)/100.2950;
T=200.59*W;
P=1.7632*pow(10,-5)*pow(W,16.0816)*ps;
}
elseif(h>86.0000&h<=91.0000)
{
W=exp((87.2848-H)/5.47);
T=186.87;
P=3.6411*pow(10,-6)*W*ps;
}
else
{
T=186.87;
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,&T);
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;
}
elsereturn0.984807753012;//=cos(10°)
}
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);
if(t>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_phi-sin_phi)+w_ie*w_ie*r*sin_psi*sin_phi*cos_phi/cos_gamma)/V;
//yDerivative[7]=1;//时间,便与调试
//yDerivative[8]=yDerivative[3];
//这一部分用于保存计算过程中的中间变量
yDerivative[7]=yDerivative[3]/g;//切向过载
yDerivative[8]=qS;//动压
yDerivative[9]=D;//阻力
yDerivative[10]=acos(cs);//sigma
yDerivative[11]=alpha;//攻角°
}
floatabove10km(double*y,floatdt){
//终止条件
//函数指针仍要调用这个函数
doubleVy=y[3]*sin(y[4]);
//if(y[1]>0.253355350897)return0.001;
//if(abs(y[4])>1.5&abs(y[5])>1.5)return0.001;
if((Vy*dt+y[0])<10000.0)
return(10000.0-y[0])/Vy;
else
return10000;
}
intreadStdTrjFile(double*stdV,double*stdGamma,double*stdnx,doub
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 航天 飞行 动力学 课程设计 飞船 再入 质点 弹道 数值 计算