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

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

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

$$ \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,y$ 为运载器在发射坐标系下位置坐标、$\phi$ 为俯仰角、$r$ 为地心距、$h$ 为飞行高度、$P$ 为总的发动机推力大小、$\alpha$ 为攻角、$g$ 为引力加速度大小、$L$ 为气动升力大小、$D$ 为气动阻力大小、$R_e$ 为地球半径、$m$ 为运载器当前质量、$m_0$ 为运载器起飞质量、$m$ 为质量秒耗量、$t$ 为飞行时间。

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

$$ 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)$ 所示。为地球引力常数。

    $$ g=\frac{\mu}{r^2}\tag{3} $$

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

$$ \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)$ 式所示。

$$ \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=\sqrt{40/(\frac{P_{0总}}{m_0*g_0}-1)}$, $P_{0总}$ 为总的地面推力,地面引力加速度 $g_0=9.8 m/s^2$。$\phi_f$、$t_2$ 为给定设计参数(如下)。

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

二、大作业任务

第一问

题目

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

$$ \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)),进行主动段轨迹仿真。给出仿真结果(飞行状态如位置、速度、高度、欧拉角等随时间变化曲线),得出两种情况下主动段终点 $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)计算结果,有:

$$ \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} $$

再因为

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

因此有:

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

再将$\Theta_k$,$r_l$,$r_k$,$v_k$带入,可得:

$$ \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)计算结果,有:

$$ \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} $$

将$x_k$,$y_t$,$R_e$代入,可得:

$$ r_k=6722273.510508056 \\ r_p=6437435.844933111 \\ r_l=6428100 \\ \quad\\ r_k\geq r_p\geq r_l \quad成立 \qquad(1) $$

再将$\Theta_k$,$r_l$,$r_k$,$v_k$带入,可得:

$$ \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)可得主动段终点的参数:

$$ \begin{aligned} & \Theta_k=23.51929491334418 \\ & \nu_k=0.5412241844023746 \\ & r_k=6642635.272693954 \\ \end{aligned} $$

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

$$ \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} $$

因此有:

$$ A\tan^2\frac{\beta_c}{2}-B\tan\frac{\beta_c}{2}+C=0 $$

再因为$\beta\geq0$,因此有:

$$ \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} $$

同理,对自由段,有:

$$ \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} $$

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

$$ \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)

解:由题意可得

导弹的最小负加速度:

$$ \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} $$

导弹的最大驻点热流:

$$ \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} $$

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

$$ \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} $$

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

(5)

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

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

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

并求得其矢积:

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

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

$$ \begin{aligned} & h_x=h\sin i\sin \Omega \\ & h_y=-h\sin i\cos \Omega \\ & h_z=h\cos i \end{aligned} $$

因此易得:

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

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

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

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

$$ \begin{aligned} & r=\frac{p}{1+e\cos(\theta-\omega)} \\ & p=h^2/\mu \end{aligned} $$

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

$$ \omega= $$

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

$$ a\cos E=ae+r\cos f \\ b\sin E=r\sin f $$

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

$$ t-t_p=\sqrt{\frac{a^3}{\mu}}(E-e\sin E) $$

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

$$ t_p=-\sqrt{\frac{a^3}{\mu}}(E-e\sin E)= $$

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

References

忘了(

最后修改:2024 年 09 月 14 日
如果觉得我的文章对你有用,请随意赞赏