侧边栏壁纸
  • 累计撰写 26 篇文章
  • 累计创建 19 个标签
  • 累计收到 4 条评论

目 录CONTENT

文章目录

【学习笔记】航天器在空中飞行的运动仿真

panedioic
2022-06-05 / 0 评论 / 0 点赞 / 319 阅读 / 5,793 字

本文仅为某次课程大作业,有诸多不严谨、遗漏之处,全文仅供参考!

一、大作业参考模型与数据

某运载器从发射点Ls采用垂直发射方式发射,假定整个飞行过程中不存在级间分离,运载器与其有效载荷一同入轨或一同再入。假定运载器主动段只在其射面内飞行,不考虑侧向运动。假定地球为均质圆球且地球周围没有风。假定飞行过程中地心坐标系为惯性坐标系且其x轴指向平春分点。主动段纵向质心动力学方程组(书中式3.70)简化为如式(1)所示形式。

v=Pcos(α)Dmgsinθθ=Pcos(α)+Lmvgcosθvx=vcosθy=vsinθ(1)ϕ=θ+αr=x2+(y+Re)2h=rRem=m0mt\begin{aligned} &v=\frac{P\cos(\alpha)-D}{m}-g\sin\theta \\ &\theta=\frac{P\cos(\alpha)+L}{mv}-g\frac{cos\theta}{v} \\ &x=v\cos\theta \\ &y=v\sin\theta \qquad\qquad\qquad\qquad\qquad\qquad\qquad(1)\\ &\phi=\theta+\alpha \\ &r=\sqrt{x^2+(y+R_e)^2} \\ &h=r-R_e \\ & m=m_0-mt \\ \end{aligned}

其中 α\alpha 为飞行速度大小、θ\theta 为速度倾角、x,yx,y 为运载器在发射坐标系下位置坐标、ϕ\phi 为俯仰角、rr 为地心距、hh 为飞行高度、PP 为总的发动机推力大小、α\alpha 为攻角、gg 为引力加速度大小、LL 为气动升力大小、DD 为气动阻力大小、ReR_e 为地球半径、mm 为运载器当前质量、m0m_0 为运载器起飞质量、mm 为质量秒耗量、tt 为飞行时间。

  • a)大气模型:大气密度采用指数模型,大气压强采用拟合模型如式 (2)(2) 所示,假定声速 vsv_s 不随飞行高度变化,vs=340m/sv_s=340m/s 为常值。

p(h)={p0e0.142h,hh92km其中p0=101325Pap(hpf),h>92km其中hpf=92km(2)p(h)= \begin{cases} p_0e^{-0.142h}, & h\le h\le 92km \quad 其中p_0=101325Pa\\ p(h_{pf}), & h>92km \quad 其中h_{pf}=92km \end{cases} \tag{2}

  • b)引力模型:采用平方反比模型,如式 (3)(3) 所示。为地球引力常数。

g=μr2(3)g=\frac{\mu}{r^2}\tag{3}

  • c)气动模型:采用拟合模型如下所示。(MM 为马赫数、CLC_L 为升力系数、CDC_D 为阻力系数、注意此拟合模型中攻角 α\alpha 单位为度)

CL(α,M)=0.0005225α2+0.03506α0.04857M+0.1577CD(α,M)=0.0001432α2+0.00558α0.01048M+0.2204\begin{aligned} &C_L(\alpha,M)=-0.0005225\alpha^2+0.03506\alpha-0.04857M+0.1577 \\ &C_D(\alpha,M)=0.0001432\alpha^2+0.00558\alpha-0.01048M+0.2204 \end{aligned}

  • d)主动段飞行程序角分段设计,如 (4)(4) 式所示。

ϕ(h)={π2tt1a(tt2)2+ϕf其中a=π2ϕf(t1t2)2t1<tt2ϕft>t2(4)\phi(h)= \begin{cases} \frac{\pi}{2} &t\le t_1 \\ a(t-t_2)^2+\phi_f \quad 其中a=\frac{\frac{\pi}{2}-\phi_f}{(t_1-t_2)^2} & t_1 < t \le t_2 \\ \phi_f & t>t_2 \end{cases} \tag{4}

其中,t=40/(P0m0g01)t=\sqrt{40/(\frac{P_{0总}}{m_0*g_0}-1)}, P0P_{0总} 为总的地面推力,地面引力加速度 g0=9.8m/s2g_0=9.8 m/s^2ϕf\phi_ft2t_2 为给定设计参数(如下)。

  • e)仿真计算采用数据
    地球半径 RE=6378.1kmR_E=6378.1km、飞行器参考面积(特征面积)为 20m220m^2、起飞质量 m0=200000kgm_0=200000kg 、飞行器空重即结构质量为 90000kg90000kg。主发动机共 33 台,单台发动机地面推力为为 1500000N1500000N、发动机喷口截面积为 0.05m20.05 m^2。主动段飞行采用耗尽关机方式(即将所有燃料全部用完发动机才关机)。
    发射点Ls经纬度分别为 λ0=60deg,ϕ0=45deg\lambda_0=60\deg, \phi_0=45\deg,地心方位角 α0=30deg\alpha_0=30\deg
    假定在发射瞬时 t=0t=0 时刻,飞行速度 v0=0.1m/sv_0=0.1m/s、速度倾角 θ0=90deg\theta_0=90\deg、运载器在发射坐标系下位置坐标 x0=0m,y0=0mx_0=0m, y_0=0m

二、大作业任务

第一问

题目

按照上述模型与数据,分别在:

Case1:m=500kg/st2=200sϕf=20deg(换成弧度)Case2:m=300kg/st2=230sϕf=5deg(换成弧度)\begin{aligned} &Case1: \quad m=500kg/s,t2=200s,\phi_f=20\deg(换成弧度) \\ &Case2: \quad m=300kg/s,t2=230s,\phi_f=5\deg(换成弧度) \end{aligned}

两种情况下通过编写主动段飞行仿真程序(数值积分式(1)),进行主动段轨迹仿真。给出仿真结果(飞行状态如位置、速度、高度、欧拉角等随时间变化曲线),得出两种情况下主动段终点 kk 处在发射坐标系下的飞行速度、速度倾角、位置坐标。

解答

由题意,弹道的主动段飞行轨迹式一阶变系数方程组,可以用数值的方法求得数值解。本题可以通过欧拉法的方式进行求解。
求解可以使用随后列出的仿真程序。在程序中输入所需的,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(νk2)cos2Θk=0.8664258403313946p=rk2vk2cos2Θk=rkνkcos2Θk=12014996\begin{aligned} & e=\sqrt{1+\nu_k(\nu_k-2)\cos^2\Theta_k}=0.8664258403313946 \\ & p=r^2_kv^2_k\cos^2\Theta_k = r_k\nu_kcos^2\Theta_k=12014996 \end{aligned}

再因为

rk=6722273.510508056mrp=6437435.844933111mrl=6428100m\begin{aligned} & r_k=6722273.510508056\, m\\ & r_p=6437435.844933111\, m \\ & r_l=6428100\, m \\ \end{aligned}

因此有:

rkrprl成立(1)r_k\geq r_p\geq r_l \quad成立 \qquad(1)

再将Θk\Theta_krlr_lrkr_kvkv_k带入,可得:

cosΘk=0.916925739930821rlrk1+2μvk2(1rl1rk)=1.0256422532299725vk2=32476863.631155732μrl1+rk/rl=59021397.478749804cosΘkrlrk1+2μvk2(1rl1rk)不成立(2)vk22μrl1+rk/rl不成立(3)\begin{aligned} & \cos\Theta_k=0.916925739930821\\ & \frac{r_l}{r_k}\sqrt{1+\frac{2\mu}{v^2_k}(\frac{1}{r_l}-\frac{1}{r_k})}=1.0256422532299725 \\ & v^2_k=32476863.63115573 \\ & \frac{2\mu r_l}{1+r_k/r_l}=59021397.478749804 \\ \\ & \cos\Theta_k \geq \frac{r_l}{r_k}\sqrt{1+\frac{2\mu}{v^2_k}(\frac{1}{r_l}-\frac{1}{r_k})}\quad不成立\qquad(2) \\ & v^2_k\geq \frac{2\mu r_l}{1+r_k/r_l}\quad不成立\qquad(3) \end{aligned}

因为(2), (3)不成立,因此此情况下能成为导弹。

case2:

根据(1)计算结果,有:

e=1+νk(νk2)cos2Θk=0.8664258403313946p=rk2vk2cos2Θk=rkνkcos2Θk=12014996.606458724\begin{aligned} & e=\sqrt{1+\nu_k(\nu_k-2)\cos^2\Theta_k}=0.8664258403313946 \\ & p=r^2_kv^2_k\cos^2\Theta_k = r_k\nu_kcos^2\Theta_k=12014996.606458724 \end{aligned}

xkx_kyty_tReR_e代入,可得:

rk=6722273.510508056rp=6437435.844933111rl=6428100rkrprl成立(1)r_k=6722273.510508056 \\ r_p=6437435.844933111 \\ r_l=6428100 \\ \quad\\ r_k\geq r_p\geq r_l \quad成立 \qquad(1)

再将Θk\Theta_krlr_lrkr_kvkv_k带入,可得:

cosΘk=0.9773518259937374rlrk1+2μvk2(1rl1rk)=0.9763775707871245vk2=110877750.456931522μrl1+rk/rl=57429289.15560445cosΘkrlrk1+2μvk2(1rl1rk)成立(2)vk22μrl1+rk/rl成立(3)\begin{aligned} & \cos\Theta_k=0.9773518259937374\\ & \frac{r_l}{r_k}\sqrt{1+\frac{2\mu}{v^2_k}(\frac{1}{r_l}-\frac{1}{r_k})}=0.9763775707871245 \\ & v^2_k=110877750.45693152 \\ & \frac{2\mu r_l}{1+r_k/r_l}=57429289.15560445 \\ \\ & \cos\Theta_k \geq \frac{r_l}{r_k}\sqrt{1+\frac{2\mu}{v^2_k}(\frac{1}{r_l}-\frac{1}{r_k})}\quad成立\qquad(2) \\ & v^2_k\geq \frac{2\mu r_l}{1+r_k/r_l}\quad成立\qquad(3) \end{aligned}

因为(1), (2), (3)均成立,因此此情况下能成为卫星。

(3).

解:由(2)可知,case1的情况可以成为导弹。此时,根据(1)可得主动段终点的参数:

Θk=23.51929491334418νk=0.5412241844023746rk=6642635.272693954\begin{aligned} & \Theta_k=23.51929491334418 \\ & \nu_k=0.5412241844023746 \\ & r_k=6642635.272693954 \\ \end{aligned}

代入被动段射程公式,可得:

A=2R(1+tan2Θk)νk(R+rk)=8125218.006125015B=2νkRtanΘk=3004693.9006598108C=νk(Rrk)=143172.88720944524\begin{aligned} & A=2R(1+\tan^2\Theta_k)-\nu_k(R+r_k)=8125218.006125015 \\ & B=2\nu_k R\tan\Theta_k=3004693.9006598108 \\ & C=\nu _k(R-r_k)=-143172.88720944524 \end{aligned}

因此有:

Atan2βc2Btanβc2+C=0A\tan^2\frac{\beta_c}{2}-B\tan\frac{\beta_c}{2}+C=0

再因为β0\beta\geq0,因此有:

tanβC2=B+B24AD2A=0.36979822647587485LKC=RβC=Rarctan(0.7084048610762707)=4518277.044430562\begin{aligned} & \tan\frac{\beta_C}{2}=\frac{B+\sqrt{B^2-4AD}}{2A}=0.36979822647587485 \\ & \quad \\ & L_{KC}=R\beta_C=R\arctan(0.7084048610762707)=4518277.044430562 \end{aligned}

同理,对自由段,有:

A=2rk(1+tan2Θk)2rkνk=8611327.040417206B=2rkνktanΘk=3129315.264447287C=0tanβE2=BA=νksinΘkcosΘk1νkcos2ΘkLKE=RβE=Rarctan(0.6971160493800195)=4446275.874550702\begin{aligned} & A=2r_k(1+\tan^2\Theta_k)-2r_k\nu_k=8611327.040417206 \\ & B=2r_k\nu_k \tan\Theta_k=3129315.264447287 \\ & C=0 \\ & \quad \\ & \tan\frac{\beta_E}{2}=\frac{B}{A}=\frac{\nu_k\sin\Theta_k\cos\Theta_k}{1-\nu_k\cos^2\Theta_k} \\ & \quad \\ & L_{KE}=R\beta_E=R\arctan(0.6971160493800195)=4446275.874550702 \end{aligned}

综上所述被动段的射程LKC=4518277.044430562L_{KC}=4518277.044430562,自由段的射程LKE=4446275.874550702L_{KE}=4446275.874550702

tanΘKOPT=νk[2Rνk(R+rk)]2Rνk4(Rrk)=0.6229555659247278tanβCmax2=νk[Rνk2(Rrk)]2[2Rνk(R+rk)]=0.4344003120021654tanΘKEOPT=1νk=0.6773299163610192tanβEmax2=12νk1νk=0.3995277421895981\begin{aligned} & \tan\Theta_{K·OPT}=\sqrt{\frac{\nu_k[2R-\nu_k(R+r_k)]}{2R\nu_k-4(R-r_k)}}=0.6229555659247278 \\ & \tan\frac{\beta_{Cmax}}{2}=\sqrt{\frac{\nu_k[R\nu_k-2(R-r_k)]}{2[2R-\nu_k(R+r_k)]}}=0.4344003120021654 \\ & \quad \\ & \tan\Theta_{KE·OPT}=\sqrt{1-\nu_k}=0.6773299163610192 \\ & \tan\frac{\beta_{Emax}}{2}=\frac{1}{2}\frac{\nu_k}{\sqrt{1-\nu_k}}=0.3995277421895981 \end{aligned}

所以,主动段终点k处的当地速度倾角不是最优速度倾角。

(4)

解:由题意可得

导弹的最小负加速度:

hm=1βln(CxSmρ0mβsinΘe)=9040.165079085044vm=vee12=3456.5257765031506vm˙βve22esinΘe=1571.31518785745\begin{aligned} & h_m=\frac{1}{\beta}\ln(-\frac{C_xS_m\rho_0}{m\beta \sin\Theta_e})=9040.165079085044 \\ & v_m=v_ee^{-\frac{1}{2}}=3456.5257765031506 \\ & \dot{v_m}\frac{\beta v^2_e}{2e}\sin\Theta_e=1571.31518785745 \end{aligned}

导弹的最大驻点热流:

hm2=1βln(3CxSmρ0mβsinΘe)=16853.908241731628vm2=vee16=4823.970321318141(qs)max=ksmβsinΘe3eCxSmve3=11825312565572.254\begin{aligned} & h_{m2}=\frac{1}{\beta}\ln(-\frac{3C_xS_m\rho_0}{m\beta \sin\Theta_e})=16853.908241731628 \\ & v_{m2}=v_ee^{-\frac{1}{6}}=4823.970321318141 \\ & (q_s)_{max}=k_s\sqrt{\frac{-m\beta\sin\Theta_e}{3eCxSm}}v^3_e=11825312565572.254 \end{aligned}

导弹的再入过程总吸热量:

vc=vee2BβsinΘeρe=161.3273896537413Q=mcfST4CxSM(ve2vc2)=9126797935641.537\begin{aligned} & v_c=v_ee^{\frac{2B}{\beta\sin\Theta_e}\rho_e}=161.3273896537413 \\ & Q=\frac{mc_f^{\prime}S_T}{4C_xS_M}(v^2_e-v^2_c)=9126797935641.537 \end{aligned}

分析最大驻点热流(qs)max(q_s)_{max}可以发现,若要降低最大驻点热流,应减小再入角Θe|\Theta_e|,但如果Θe|\Theta_e|过小,又会增加飞行的时间,使总吸热量增加。因此,调整时需要兼顾总吸热量。合理的选择一个再入角Θe\Theta_e,是弹道再入的一个重要问题。

(5)

解:由(2)可知,case2的情况可以成为卫星。此时,根据(1)可得主动段终点的参数:

由上述参数可求出弹道关机点时的卫星地心矩矢量r速度矢量v分别为:

r=[369797.99601624296122.0692899944279.81467744]Tv=[4513.921402073078.239390851620.43804561]T\begin{aligned} & r=[-369797.99601624 \quad 296122.06928999 \quad -44279.81467744]^T \\ & v=[-4513.92140207 \quad 3078.23939085 \quad -1620.43804561]^T \end{aligned}

并求得其矢积:

h=r×v=[3.43543597e+083.99359139e+081.98344988e+08]\begin{aligned} h&=r\times v \\ &=[-3.43543597e+08 \quad -3.99359139e+08 \quad 1.98344988e+08] \end{aligned}

由卫星轨道与赤道惯性坐标系之间的关系可得动量矩h的三个分量与iiΩ\Omega的关系:

hx=hsinisinΩhy=hsinicosΩhz=hcosi\begin{aligned} & h_x=h\sin i\sin \Omega \\ & h_y=-h\sin i\cos \Omega \\ & h_z=h\cos i \end{aligned}

因此易得:

i=arccos(hzh)=1.2106983360402104radΩ=arctan(hzhy)=0.46097059130336926rad\begin{aligned} & i=\arccos(\frac{h_z}{h})=1.2106983360402104\,rad \\ & \Omega=\arctan(\frac{-h_z}{h_y})=0.46097059130336926\,rad \end{aligned}

由开普勒定律可知,卫星绕地球的运动轨迹为一椭圆。因此易得:

e=0.8664258403313946a=h2μ(1e2)=\begin{aligned} & e=0.8664258403313946 \\ & a=\frac{h^2}{\mu(1-e^2)}= \end{aligned}

再因为卫星轨道方程的标准极坐标形式:

r=p1+ecos(θω)p=h2/μ\begin{aligned} & r=\frac{p}{1+e\cos(\theta-\omega)} \\ & p=h^2/\mu \end{aligned}

将关机点参数r=,θ=0r=, \theta=0代入,可解得近地点幅角:

ω=\omega=

由几何关系可得出卫星轨道的偏近点角EE与真近点角ff之间的关系:

acosE=ae+rcosfbsinE=rsinfa\cos E=ae+r\cos f \\ b\sin E=r\sin f

已知卫星的位置与时间关系的开普勒方程:

ttp=a3μ(EesinE)t-t_p=\sqrt{\frac{a^3}{\mu}}(E-e\sin E)

将关机点参数f=ω,t=0f=-\omega, t=0代入,可解得卫星经过近地点的时刻tpt_p

tp=a3μ(EesinE)=t_p=-\sqrt{\frac{a^3}{\mu}}(E-e\sin E)=

至此,给出了卫星运动方程的6个积分常数。

References

忘了(

0

评论区