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

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

某运载器从发射点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

忘了（