本文仅为某次课程大作业,有诸多不严谨、遗漏之处,全文仅供参考!
一、大作业参考模型与数据
某运载器从发射点Ls采用垂直发射方式发射,假定整个飞行过程中不存在级间分离,运载器与其有效载荷一同入轨或一同再入。假定运载器主动段只在其射面内飞行,不考虑侧向运动。假定地球为均质圆球且地球周围没有风。假定飞行过程中地心坐标系为惯性坐标系且其x轴指向平春分点。主动段纵向质心动力学方程组(书中式3.70)
简化为如式(1)
所示形式。
v=mPcos(α)−D−gsinθθ=mvPcos(α)+L−gvcosθx=vcosθy=vsinθ(1)ϕ=θ+αr=x2+(y+Re)2h=r−Rem=m0−mt
其中 α 为飞行速度大小、θ 为速度倾角、x,y 为运载器在发射坐标系下位置坐标、ϕ 为俯仰角、r 为地心距、h 为飞行高度、P 为总的发动机推力大小、α 为攻角、g 为引力加速度大小、L 为气动升力大小、D 为气动阻力大小、Re 为地球半径、m 为运载器当前质量、m0 为运载器起飞质量、m 为质量秒耗量、t 为飞行时间。
- a)大气模型:大气密度采用指数模型,大气压强采用拟合模型如式 (2) 所示,假定声速 vs 不随飞行高度变化,vs=340m/s 为常值。
p(h)={p0e−0.142h,p(hpf),h≤h≤92km其中p0=101325Pah>92km其中hpf=92km(2)
- b)引力模型:采用平方反比模型,如式 (3) 所示。为地球引力常数。
g=r2μ(3)
- c)气动模型:采用拟合模型如下所示。(M 为马赫数、CL 为升力系数、CD 为阻力系数、注意此拟合模型中攻角 α 单位为度)
CL(α,M)=−0.0005225α2+0.03506α−0.04857M+0.1577CD(α,M)=0.0001432α2+0.00558α−0.01048M+0.2204
- d)主动段飞行程序角分段设计,如 (4) 式所示。
ϕ(h)=⎩⎪⎪⎨⎪⎪⎧2πa(t−t2)2+ϕf其中a=(t1−t2)22π−ϕfϕft≤t1t1<t≤t2t>t2(4)
其中,t=40/(m0∗g0P0总−1), P0总 为总的地面推力,地面引力加速度 g0=9.8m/s2。ϕf、t2 为给定设计参数(如下)。
- e)仿真计算采用数据
地球半径 RE=6378.1km、飞行器参考面积(特征面积)为 20m2、起飞质量 m0=200000kg 、飞行器空重即结构质量为 90000kg。主发动机共 3 台,单台发动机地面推力为为 1500000N、发动机喷口截面积为 0.05m2。主动段飞行采用耗尽关机方式(即将所有燃料全部用完发动机才关机)。
发射点Ls经纬度分别为 λ0=60deg,ϕ0=45deg,地心方位角 α0=30deg。
假定在发射瞬时 t=0 时刻,飞行速度 v0=0.1m/s、速度倾角 θ0=90deg、运载器在发射坐标系下位置坐标 x0=0m,y0=0m。
二、大作业任务
第一问
题目
按照上述模型与数据,分别在:
Case1:m=500kg/s,t2=200s,ϕf=20deg(换成弧度)Case2:m=300kg/s,t2=230s,ϕf=5deg(换成弧度)
两种情况下通过编写主动段飞行仿真程序(数值积分式(1)),进行主动段轨迹仿真。给出仿真结果(飞行状态如位置、速度、高度、欧拉角等随时间变化曲线),得出两种情况下主动段终点 k 处在发射坐标系下的飞行速度、速度倾角、位置坐标。
解答
由题意,弹道的主动段飞行轨迹式一阶变系数方程组,可以用数值的方法求得数值解。本题可以通过欧拉法的方式进行求解。
求解可以使用随后列出的仿真程序。在程序中输入所需的,t2,的值以及合适的步长,即可求出所需的数值解。
以下为在设定不同的步长时仿真所的出的case1的x坐标值、所花费的时间与步长的关系图:
观察可发现,当步长设定为0.001时计算花费了。。的时间,且在此之后所得出的数值差距不大。因此之后的仿真均采用0.001作为步长进行计算。
经过计算,可得在case1与case2下主动段终点k处的飞行速度、飞行倾角、位置坐标分别为:
并绘出飞行状态随时间的变化曲线(对于每个case的曲线分别为:速度-时间变化曲线、飞行高度-时间变化曲线、y-x曲线、攻角-时间变化曲线、速度倾角-时间变化曲线、俯仰角-时间变化曲线):
附:计算所需的程序:
# homework.py
import numpy as np
import matplotlib.pyplot as plt
# basic constants
PI = 3.14159265
beginMass = 200000 # mkg
beginGravity = 9.8
structureMass = 90000
Sref = 20
Propulsion = 1500000 * 3 # m3 * 1.5e6 N
baseDensity = 1.125 # mkg/m^3
baseAcousitc = 340 # mm/s
Rearth = 6378100 # m
mu = 398600000000000 # m^3/s^2
def getGravity(radius):
return mu / radius / radius
def getAtmosphericDensity(height):
return baseDensity * np.exp(-1.406 * 10 ** (-4) * height)
def getAerodynamicDrag(height, alpha, velocity):
density = getAtmosphericDensity(height)
mach = velocity / baseAcousitc
alpha_deg = alpha * 180 / 3.14159265
coefficient = 0.0001432 * alpha_deg ** 2 + 0.00558 * alpha_deg - 0.01048 * mach + 0.2204
return 0.5 * coefficient * density * velocity **2 * Sref
def getAerodynamicLift(height, alpha, velocity):
density = getAtmosphericDensity(height)
mach = velocity / baseAcousitc
alpha_deg = alpha * 180 / 3.14159265
coefficient = -0.0005225 * alpha_deg ** 2 + 0.03506 * alpha_deg - 0.04857 * mach + 0.1577
return 0.5 * coefficient * density * velocity **2 * Sref
def getPhi(time, t1, t2, phiF):
if(time < t1):
return PI / 2
elif (time < t2):
a = (PI / 2 - phiF) / (t1-t2)**2
return a * (time-t2)**2 + phiF
else :
return phiF
def run(dm, t2, phiF, stepLength):
# calculate constant
t1 = np.sqrt(40/(Propulsion / (beginMass * beginGravity) - 1))
# calculate variables
nowTime = 0.0
mass = beginMass
v = 0.1
x = 0
y = 0
theta = PI / 2
dataSet = []
cnt = 0
while(True):
lengthToEarthCenter = np.sqrt(x ** 2 + (y + Rearth) ** 2)
height = lengthToEarthCenter - Rearth
phi = getPhi(nowTime, t1, t2, phiF)
alpha = phi - theta
dv = (Propulsion * np.cos(alpha) - getAerodynamicDrag(height, alpha, v)) / mass - getGravity(lengthToEarthCenter) * np.sin(theta)
dtheta = (Propulsion * np.sin(alpha) + getAerodynamicLift(height, alpha, v)) / (mass * v) - getGravity(lengthToEarthCenter) * np.cos(theta) / v
v += dv * stepLength
theta += dtheta * stepLength
x += v * np.cos(theta) * stepLength
y += v * np.sin(theta) * stepLength
mass -= dm * stepLength
nowTime += stepLength
cnt += 1
if(mass <= structureMass):
break
if(cnt % 1000 == 0):
dataSet.append([nowTime, alpha, theta, phi, x, y, height, v])
print('计算结束!')
print('位置(xk, yk, zk) = ('+str(x/1000)+'km, '+str(y/1000)+'km, 0km)、速度vk='+str(v)+'、速度倾角θk='+str((theta+np.arctan(x/(y+Rearth)))*180/3.14159265)+'°')
return dataSet
def showData(datas):
for data in datas:
time = [el[0] for el in data]
alpha = [el[1] for el in data]
theta = [el[2] for el in data]
phi = [el[3] for el in data]
x = [el[4] for el in data]
y = [el[5] for el in data]
height = [el[6] for el in data]
v = [el[7] for el in data]
# 使用子图
plt.figure()
plt.suptitle('The plot of case '+str(cnt), size=16) # 整张图的总标题
plt.subplot(2, 3, 1) # (行,列,活跃区)
plt.title('velocity - time') # 子图的小标题
plt.plot(time, v, ls="-", lw=2)
plt.subplot(2, 3, 2)
plt.title('height - time') # 子图的小标题
plt.plot(time, height, ls="-", lw=2)
plt.subplot(2, 3, 3)
plt.title('y - x') # 子图的小标题
plt.plot(x, y, ls="-", lw=2)
plt.subplot(2, 3, 4)
plt.title('alpha - time') # 子图的小标题
plt.plot(time, alpha, ls="-", lw=2)
plt.subplot(2, 3, 5)
plt.title('theta - time') # 子图的小标题
plt.plot(time, theta, ls="-", lw=2)
plt.subplot(2, 3, 6)
plt.title('phi - time') # 子图的小标题
plt.plot(time, phi, ls="-", lw=2)
plt.show()
print('case1:')
case1 = run(500, 200, 20.0 / 180.0 * 3.14159265, 10**(0))
print('case2:')
case2 = run(300, 230, 5.0 / 180.0 * 3.14159265, 0.35)
showData([case1, case2])
程序的输出为:
case1:
位置(xk, yk, zk) = (403.4326437092588km, 252.27295090009528km, 0km)、速度vk=5698.847570444022、速度倾角θk=23.51929491334418°
case2:
位置(xk, yk, zk) = (1449.6315824637384km, 228.72349043694848km, 0km)、速度vk=10529.850447985076、速度倾角θk=12.21735847248481°
第二问
case1:
根据(1)计算结果,有:
e=1+νk(νk−2)cos2Θk=0.8664258403313946p=rk2vk2cos2Θk=rkνkcos2Θk=12014996
再因为
rk=6722273.510508056mrp=6437435.844933111mrl=6428100m
因此有:
rk≥rp≥rl成立(1)
再将Θk,rl,rk,vk带入,可得:
cosΘk=0.916925739930821rkrl1+vk22μ(rl1−rk1)=1.0256422532299725vk2=32476863.631155731+rk/rl2μrl=59021397.478749804cosΘk≥rkrl1+vk22μ(rl1−rk1)不成立(2)vk2≥1+rk/rl2μrl不成立(3)
因为(2), (3)不成立,因此此情况下能成为导弹。
case2:
根据(1)计算结果,有:
e=1+νk(νk−2)cos2Θk=0.8664258403313946p=rk2vk2cos2Θk=rkνkcos2Θk=12014996.606458724
将xk,yt,Re代入,可得:
rk=6722273.510508056rp=6437435.844933111rl=6428100rk≥rp≥rl成立(1)
再将Θk,rl,rk,vk带入,可得:
cosΘk=0.9773518259937374rkrl1+vk22μ(rl1−rk1)=0.9763775707871245vk2=110877750.456931521+rk/rl2μrl=57429289.15560445cosΘk≥rkrl1+vk22μ(rl1−rk1)成立(2)vk2≥1+rk/rl2μrl成立(3)
因为(1), (2), (3)均成立,因此此情况下能成为卫星。
(3).
解:由(2)可知,case1的情况可以成为导弹。此时,根据(1)可得主动段终点的参数:
Θk=23.51929491334418νk=0.5412241844023746rk=6642635.272693954
代入被动段射程公式,可得:
A=2R(1+tan2Θk)−νk(R+rk)=8125218.006125015B=2νkRtanΘk=3004693.9006598108C=νk(R−rk)=−143172.88720944524
因此有:
Atan22βc−Btan2βc+C=0
再因为β≥0,因此有:
tan2βC=2AB+B2−4AD=0.36979822647587485LKC=RβC=Rarctan(0.7084048610762707)=4518277.044430562
同理,对自由段,有:
A=2rk(1+tan2Θk)−2rkνk=8611327.040417206B=2rkνktanΘk=3129315.264447287C=0tan2βE=AB=1−νkcos2ΘkνksinΘkcosΘkLKE=RβE=Rarctan(0.6971160493800195)=4446275.874550702
综上所述被动段的射程LKC=4518277.044430562,自由段的射程LKE=4446275.874550702。
tanΘK⋅OPT=2Rνk−4(R−rk)νk[2R−νk(R+rk)]=0.6229555659247278tan2βCmax=2[2R−νk(R+rk)]νk[Rνk−2(R−rk)]=0.4344003120021654tanΘKE⋅OPT=1−νk=0.6773299163610192tan2βEmax=211−νkνk=0.3995277421895981
所以,主动段终点k处的当地速度倾角不是最优速度倾角。
(4)
解:由题意可得
导弹的最小负加速度:
hm=β1ln(−mβsinΘeCxSmρ0)=9040.165079085044vm=vee−21=3456.5257765031506vm˙2eβve2sinΘe=1571.31518785745
导弹的最大驻点热流:
hm2=β1ln(−mβsinΘe3CxSmρ0)=16853.908241731628vm2=vee−61=4823.970321318141(qs)max=ks3eCxSm−mβsinΘeve3=11825312565572.254
导弹的再入过程总吸热量:
vc=veeβsinΘe2Bρe=161.3273896537413Q=4CxSMmcf′ST(ve2−vc2)=9126797935641.537
分析最大驻点热流(qs)max可以发现,若要降低最大驻点热流,应减小再入角∣Θe∣,但如果∣Θe∣过小,又会增加飞行的时间,使总吸热量增加。因此,调整时需要兼顾总吸热量。合理的选择一个再入角Θe,是弹道再入的一个重要问题。
(5)
解:由(2)可知,case2的情况可以成为卫星。此时,根据(1)可得主动段终点的参数:
由上述参数可求出弹道关机点时的卫星地心矩矢量r速度矢量v分别为:
r=[−369797.99601624296122.06928999−44279.81467744]Tv=[−4513.921402073078.23939085−1620.43804561]T
并求得其矢积:
h=r×v=[−3.43543597e+08−3.99359139e+081.98344988e+08]
由卫星轨道与赤道惯性坐标系之间的关系可得动量矩h的三个分量与i,Ω的关系:
hx=hsinisinΩhy=−hsinicosΩhz=hcosi
因此易得:
i=arccos(hhz)=1.2106983360402104radΩ=arctan(hy−hz)=0.46097059130336926rad
由开普勒定律可知,卫星绕地球的运动轨迹为一椭圆。因此易得:
e=0.8664258403313946a=μ(1−e2)h2=
再因为卫星轨道方程的标准极坐标形式:
r=1+ecos(θ−ω)pp=h2/μ
将关机点参数r=,θ=0代入,可解得近地点幅角:
ω=
由几何关系可得出卫星轨道的偏近点角E与真近点角f之间的关系:
acosE=ae+rcosfbsinE=rsinf
已知卫星的位置与时间关系的开普勒方程:
t−tp=μa3(E−esinE)
将关机点参数f=−ω,t=0代入,可解得卫星经过近地点的时刻tp:
tp=−μa3(E−esinE)=
至此,给出了卫星运动方程的6个积分常数。
References
忘了(
评论区