一 前言

飞行力学（flight mechanics）是研究空间中（尤其是大气内的）飞行器运动状态的学科。对了解飞行器飞行原理、进行飞行控制等诸多方面有着至关重要的作用。

目前，国内的爱好者制作的模型飞机的翼展通常在1.5m2以下，属于小型机型，机体较为灵活，航速相对较小。对其进行飞行力学分析较为容易。

对于固定翼航模机领域，国外的技术相对较为成熟。例如，2017年左右来自美国的爱好者Spencer Lisenby制作的滑翔机，在山顶背风面的两股不同流速的气流下无动力飞行。其航速可被加速到惊人的245m/s。对此我坛的大佬三水合番也发过相应的帖子（文献[1]）动态滑翔—— 245m/s 的无动力滑翔机。

本文将通过飞行力学的角度对大气内飞行器进行运动分析和理论力学建模，同时展开飞行过程中机体姿态控制的理论分析。

（图-1 Spencer Lisenby发布在YouTube上的飞机）

二 无控固定翼飞行器飞行的动力学描述

对于在大气内飞行的物体，要想全面描述其运动是件非常复杂的事。所以在本文中我们将以在空间中无侧向来风的这种情况作为理想条件，对飞行器进行力学分析。

由于飞行器的姿态会影响飞行的气动特征，故不能视飞行器为质点。在此将其看作一般刚体的自由运动并对刚体的六个自由度进行分析。由于无侧向来风且飞行器为无控飞行状态，故只考虑两个平动度和一个转动度。采用Newton力学体系下一般刚体自由运动的两条本构方程，

$m\ddot{r} = \sum{F}$ (1)

$\frac{\mathrm{d} J}{\mathrm{d} t} = \sum{M}$ (2)

做机体运动受力图（图-2），将过飞行器质心的水平连线作为飞行基线，过质心且和机体轴线共线的直线作为飞行视线，记两直线的夹角$\beta$为飞行俯仰角。

以地面为惯性系，建立空间直角坐标系，y轴垂直于地面向上。由理想条件略去z轴方向运动。

（图-2 机体受力图）

为了让飞行更加的稳定，应先确定飞行器质心和压心的位置。如图-3，随着起飞后飞行速度的增大，升力L的力矩也随之增大。为稳定机身，要将质心适当的前移。以让Ft、D等气动作用力产生的恢复力矩调整姿态并保持平稳。所以，质心应保证位于压心前30%到40%相对机体长度。如果质心过于靠后，机体稳定性将下降，但若质心过于靠前，机体灵活性又将下降，不易被控制。

（图-3 机体静力分析图）

那么，现在对机体先进行平动分析，由图-2中的受力列得刚体平动方程，

$m\frac{\mathrm{d} v_{y}}{\mathrm{d} t} = -F_{t}\sin \beta -D\cos \alpha -L_{0}\sin \beta$ (3)

$m\frac{\mathrm{d} v_{x}}{\mathrm{d} t} = F_{t}\cos \beta -D\sin \alpha -L_{0}\cos \beta - mg$ (4)

式中$\alpha$角为轨道切线偏转角，通常在弹道系下和$\beta$等价。下面对此时的阻力进行分析。考虑飞行雷诺数，

$Re = \frac{\rho v l}{\mu }$

由实验经验我们知道当流体做层流流动时（Re<2000），其所受阻力满足斯托克斯流体阻力定理,即其所受气体阻力和来流表观流速的一次方成正比。然而，此时飞行器的飞行雷诺数通常都大于甚至远大于该值，故此时采用湍流阻力的半经验公式,

$D=\frac{1}{2}C_{D}\rho A_{e}v^2$ (5)

在该理想模型下，机体的姿态将会直接影响到阻力系数的值。故借鉴文献[2]的做法，在此引入气动导数$C_{D}^\beta$、$C_{L}^\beta$。设其满足函数关系$C_{D}^\beta (\beta)$。将$C_{D}^\beta (\beta)$ 写作幂级数形式，

对（5）式稍加变形可得，

$C_{D}^\beta = \frac{2D}{\beta^2\rho A_{e} v^2} - \frac{C_{D0}}{\beta^2}$ (6)

采用风洞试验或有限元方法测量阻力D，代入（6）式，得到满足$C_{\beta} (\beta)$的数据点。将每个数值点的横纵坐标画为散点图，通过构造线性方程组来回归求系数an的值，

$\begin{bmatrix} &1 &\beta_{1} &\beta_{1}^2 &... &\beta_{1}^n\\ &1 &\beta_{2} &\beta_{2}^2 &... &\beta_{2}^n\\ & & &\vdots & &\\ &1 &\beta_{n} &\beta_{n}^2 &...&\beta_{n}^n \end{bmatrix} \begin{bmatrix} a_{0}\\ a_{1}\\ \vdots\\a_{n} \end{bmatrix} = \begin{bmatrix} C_{D1}^\beta\\ C_{D2}^\beta\\ \vdots\\C_{Dn}^\beta \end{bmatrix}$ (7)

用Gaussian消元法或Cramer法则原理编制的计算机程序解算（7）式，将系数an回代入（6）式即可。一般对于实际模型求解，取二阶精度即可。对于（3）（4）式中的升力L，满足zhukovsky条件，

$L = \rho v \Gamma$

令$C_{L}^\beta = \frac{\rho \Gamma }{\beta } -\frac{C_{L0}}{\beta}$,可以用上述相同的方法得到$C_{L}^\beta$。同理也可解决$F_{t}$项。因此，机体的质心平动方程可化为，

$m\frac{\mathrm{d} v_{x}}{\mathrm{d} t} = -\frac{1}{2}\rho v^2 A_{t}(C_{t0}+C_{t}^\beta \beta^2)\sin \beta - \frac{1}{2}\rho v^2 A_{e}(C_{D0}+C_{D}^\beta \beta^2)\cos \beta - \frac{1}{2}\rho v^2 A_{e}(C_{L0}+C_{L}^\beta \beta)\sin \beta$ (8)

$m\frac{\mathrm{d} v_{y}}{\mathrm{d} t} = \frac{1}{2}\rho v^2 A_{t}(C_{t0}+C_{t}^\beta \beta^2)\cos \beta - \frac{1}{2}\rho v^2 A_{e}(C_{D0}+C_{D}^\beta \beta^2)\sin \beta + \frac{1}{2}\rho v^2 A_{e}(C_{L0}+C_{L}^\beta \beta)\cos \beta - mg$ (9)

接下来分析机体的转动，同样由图-2，得到合力矩的表达，

$\sum M = F_{t}\times r_{t} +L \times r_{L} - D\times r_{D}$

对角动量定理（2）式稍加变形，

$\sum M = I\ddot{\beta } = I\ddot{\omega_{\beta } }$ (10)

式中I为飞行器的转动惯量，假设力总垂直于机身轴线。联立（8）（9）（10）式，得飞行器在无侧向来风和不考虑地球科里奥利效应的无动力飞行刚体动力学方程组，

$\left\{ \begin{aligned} m\frac{\mathrm{d} v_{x}}{\mathrm{d} t} &= -\frac{1}{2}\rho v^2 A_{t}(C_{t0}+C_{t}^\beta \beta^2)\sin \beta - \frac{1}{2}\rho v^2 A_{e}(C_{D0}+C_{D}^\beta \beta^2)\cos \beta - \frac{1}{2}\rho v^2 A_{e}(C_{L0}+C_{L}^\beta \beta)\sin \beta \\ m\frac{\mathrm{d} v_{y}}{\mathrm{d} t} &= \frac{1}{2}\rho v^2 A_{t}(C_{t0}+C_{t}^\beta \beta^2)\cos \beta - \frac{1}{2}\rho v^2 A_{e}(C_{D0}+C_{D}^\beta \beta^2)\sin \beta + \frac{1}{2}\rho v^2 A_{e}(C_{L0}+C_{L}^\beta \beta)\cos \beta - mg \\ I\frac{\mathrm{d} y}{\mathrm{d} x} &= \frac{1}{2} \rho v^2 A_{t} r_{t}(m_{t}^\beta \beta + m_{t}^\omega \omega ) + \frac{1}{2} \rho v^2 A_{e} r_{L}(m_{L}^\beta \beta + m_{L}^\omega \omega ) - \frac{1}{2} \rho v^2 A_{e} r_{D}(m_{D}^\beta \beta + m_{D}^\omega \omega ) \end{aligned} \right.$ (11)

对于该式中，$v = \sqrt{v_{x}^2+v_{y}^2}$、假定力矩气动导数和飞行俯仰角和转动角速度均相关。在大气中，空气的密度$\rho$将随高度h和温度T的变化而变化。但由于本文研究的对象为小型固定翼飞行器，其高度的影响对飞行的影响并不大。故假设密度不随高度变化，即$\rho(h,T) = const$。对（11）式可采用计算机运行4RK算法求解，再用数值积分法求出其轨道方程。

然而，在实际情况下飞行器为了达到更远的航程，往往会带上助推器进行动力飞行。那么整个飞行过程将分为动力飞行段和惯性滑翔段两个飞行阶段。同本节开头的质、压心稳定性分析一样。对于带推力的飞行器，由图-4所示。为了使飞行稳定，应让质心压心和推力中心三点共线于机身轴线上。且发动机轴线应平行于机身轴线，其次多台发动机应沿机身轴线对称安装。

(图-4 三心分布)

对于带动力的飞行器其动力学方程写为以下（12）式，

$\left\{ \begin{aligned} m\frac{\mathrm{d} v_{x}}{\mathrm{d} t} &= T(t)\cos \beta -\frac{1}{2}\rho v^2 A_{t}(C_{t0}+C_{t}^\beta \beta^2)\sin \beta - \frac{1}{2}\rho v^2 A_{e}(C_{D0}+C_{D}^\beta \beta^2)\cos \beta - \frac{1}{2}\rho v^2 A_{e}(C_{L0}+C_{L}^\beta \beta)\sin \beta \\ m\frac{\mathrm{d} v_{y}}{\mathrm{d} t} &= T(t)\sin \beta + \frac{1}{2}\rho v^2 A_{t}(C_{t0}+C_{t}^\beta \beta^2)\cos \beta - \frac{1}{2}\rho v^2 A_{e}(C_{D0}+C_{D}^\beta \beta^2)\sin \beta + \frac{1}{2}\rho v^2 A_{e}(C_{L0}+C_{L}^\beta \beta)\cos \beta - mg \\ I\frac{\mathrm{d} y}{\mathrm{d} x} &= \frac{1}{2} \rho v^2 A_{t} r_{t}(m_{t}^\beta \beta + m_{t}^\omega \omega ) + \frac{1}{2} \rho v^2 A_{e} r_{L}(m_{L}^\beta \beta + m_{L}^\omega \omega ) - \frac{1}{2} \rho v^2 A_{e} r_{D}(m_{D}^\beta \beta + m_{D}^\omega \omega ) \end{aligned} \right.$ (12)

当然，采用传统的Newton力学体系免不了繁琐的矢量计算。为了简化数学模型、避免理想约束内力出现和方便程式化。也可采用Language力学进行建模，

$\frac{\mathrm{d}}{\mathrm{d} t}(\frac{\partial L}{\partial \dot{q_{j}}} )-\frac{\partial L}{\partial q_{j}} = Q_{j}$ $j = 1,2,3,...$ (13)

取上面的三个独立坐标作为广义坐标，

$\left\{\begin{matrix} q_{1} = x\\ q_{2} = y\\q_{3} = \beta \end{matrix}\right.$

先明确机体上非保守力在弹道系B上的状态空间表达，

$R = \begin{bmatrix} -D\\L+F_{t} \end{bmatrix}_{B}$ (14)

沿地面惯性系A变换，

$R = \begin{bmatrix}- D \cos \beta - L \sin \beta - F_{t} \sin \beta \\- D \sin \beta +L \cos \beta + F_{t} \cos \beta \end{bmatrix}_{A}$ (15)

其转动合力矩同上述Newton力学解法中的一样，

$\sum M_{R} = F_{t} \times r_{t} + L \times r_{L} - D \times r_{D}$ (16)

得到各非保守力广义力表达集，

$\left\{ \begin{aligned}Q_{1} = - D \cos \beta - L \sin \beta - F_{t} \sin \beta \\Q_{2} = - D \sin \beta +L \cos \beta + F_{t} \cos \beta \\Q_{3} = F_{t} \times r_{t} + L \times r_{L} - D \times r_{D} \end{aligned} \right.$ (17)

再分析保守力的能量关系，由上述选取的广义坐标，不难得到机体平动和转动的动能，

$T_{1} = \frac{1}{2} m (\dot x^2 + \dot y^2)$ (18)

$T_{2} = \frac{1}{2} I \dot{\beta ^2}$ (19)

对质心，其重力势能为，

$V = mgy$ (20)

由（18）（19）（20）式，可得拉格朗日函数L，

$L = T_{1} + T_{2} - V = \frac{1}{2} m (\dot x^2 + \dot y^2) + \frac{1}{2} I \dot{\beta ^2} - mgy$ (21)

由之前Newton力学解法中的条件，联立（13）（17）（21）求解该Euler-lagrange方程，稍加化简可得，

$\begin{cases} m\ddot{x} = -D\cos\beta - L\sin\beta - F_t\sin\beta \\ m\ddot{y} + mg = -D\sin\beta + L\cos\beta + F_t\cos\beta \\ I\ddot{\beta} = F_tr_t + Lr_L - Dr_D \end{cases}$ (22)

代入气动参数、推力函数和方程（5）（6）式，可以得到和（12）式一样的结果。

利用python中的Scipy、matplotlib等库，使用Fourth-order Runge Kutta算法（4RK）对（12）式进行求解，得到各坐标和时间的关系。可以再对其进行Simpson积分法进行数值积分，得到飞行器飞行轨迹图。

在此以简单的线性条件为初值条件，以下是笔者给解算程序设置的初值条件：

m = 1.0                # 质量 [kg]
I = 0.05               # 转动惯量 [kg·m²]
A_t = 0.02             # 推力面面积 [m²]
A_e = 0.1              # 升/阻力面面积 [m²]
r_t = 0.05             # 推力臂长 [m]
r_L = 0.1              # 升力臂长 [m]
r_D = 0.1              # 阻力臂长 [m]
rho = 1.225            # 空气密度 [kg/m³]

# 气动导数
Ct0 = 0.01
Ct_beta = 0.1
CD0 = 0.05
CD_beta = 0.05
CL0 = 0.2
CL_beta = 0.1
mt_beta = 0.01
mt_omega = 0.005
mL_beta = 0.01
mL_omega = 0.005
mD_beta = 0.005
mD_omega = 0.002

#推力函数（假定，python语法）Net
T0 = 5.0   # 初始推力 [N]
def thrust(t):
    return T0

# 初始条件
Y0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  # [vx, vy, omega, beta, x, y]

具体求解算法源代码如下，

import numpy as np                                    #导库
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 运动参数定义
m = 1.0
I = 0.05
A_t = 0.02
A_e = 0.1
r_t = 0.05
r_L = 0.1
r_D = 0.1
rho = 1.225
g = 9.81

# 气动导数定义
Ct0 = 0.01
Ct_beta = 0.1
CD0 = 0.05
CD_beta = 0.05
CL0 = 0.2
CL_beta = 0.1
mt_beta = 0.01
mt_omega = 0.005
mL_beta = 0.01
mL_omega = 0.005
mD_beta = 0.005
mD_omega = 0.002

# 推力函数（恒定）
T0 = 5.0
def thrust(t):
    return T0

# 定义ODE方程函数
def dydt(t, Y):
    vx, vy, omega, beta, x, y = Y
    v = np.sqrt(vx**2 + vy**2)
    if v == 0:
        v = 1e-6  # 除0错误排除

    # 推力函数项
    T = thrust(t)

    # 力项
    F_tx = 0.5 * rho * v**2 * A_t * (Ct0 + Ct_beta * beta**2) * np.sin(beta)
    F_dx = 0.5 * rho * v**2 * A_e * (CD0 + CD_beta * beta**2) * np.cos(beta)
    F_lx = 0.5 * rho * v**2 * A_e * (CL0 + CL_beta * beta) * np.sin(beta)
    F_ty = 0.5 * rho * v**2 * A_t * (Ct0 + Ct_beta * beta**2) * np.cos(beta)
    F_dy = 0.5 * rho * v**2 * A_e * (CD0 + CD_beta * beta**2) * np.sin(beta)
    F_ly = 0.5 * rho * v**2 * A_e * (CL0 + CL_beta * beta) * np.cos(beta)

    dvx = (T * np.cos(beta) - F_tx - F_dx - F_lx) / m
    dvy = (T * np.sin(beta) + F_ty - F_dy + F_ly - m * g) / m

    # 力矩项
    M_t = 0.5 * rho * v**2 * A_t * r_t * (mt_beta * beta + mt_omega * omega)
    M_L = 0.5 * rho * v**2 * A_e * r_L * (mL_beta * beta + mL_omega * omega)
    M_D = 0.5 * rho * v**2 * A_e * r_D * (mD_beta * beta + mD_omega * omega)
    domega = (M_t + M_L - M_D) / I

    # 角度变化率（假设 dbeta/dt = omega）
    dbeta = omega

    # 位移导数
    dx = vx
    dy = vy

    return [dvx, dvy, domega, dbeta, dx, dy]

# 初始条件
Y0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  # [vx, vy, omega, beta, x, y]

# 时间范围
t_span = [0, 10]
t_eval = np.linspace(*t_span, 1000)

# 解算
sol = solve_ivp(dydt, t_span, Y0, t_eval=t_eval, method='RK45')

# 结果
t = sol.t
vx, vy, omega, beta, x, y = sol.y

# 绘图
plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)
plt.plot(t, vx, label=r'$v_x$')
plt.xlabel('Time [s]')
plt.ylabel(r'$v_x$ [m/s]')
plt.grid()

plt.subplot(2, 2, 2)
plt.plot(t, vy, 'g', label=r'$v_y$')
plt.xlabel('Time [s]')
plt.ylabel(r'$v_y$ [m/s]')
plt.grid()

plt.subplot(2, 2, 3)
plt.plot(t, omega, 'r', label=r'$\omega$')
plt.xlabel('Time [s]')
plt.ylabel(r'$\omega$ [rad/s]')
plt.grid()

plt.subplot(2, 2, 4)
plt.plot(x, y)
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.title('Trajectory')
plt.grid()

plt.tight_layout()
plt.show()

经过运行后，可以得到运算结果如下图。

（图-5 笔者的运行结果，前三个为q-t图像，最后一个为运动轨迹图）

由于实际的飞行环境十分复杂，往往需要考虑到更多的不确定因素，如变质量、变推力和气流扰动。故本文在此只进行理想条件的近似计算。

三 姿态航向调节原理

对于固定翼飞行器姿态的调整和控制，目前可供选择的方法有很多。本文以控制面(Control surface)的方法进行整体调整为例，将详细的进一步解释该原理的作用。

如图-6，当控制面沿一轴旋转后，迎风部和气流作用产生侧向偏向力从而改变姿态。对于固定翼飞机而言，姿态的改变往往伴随着航向的变化。

（图-6 控制面控制原理简示图）

对于产生的侧向力的描述，有着类似于（5）式的表达形式，

$Y = \frac{1}{2} C_{Y}\rho A_{t} v^2$ (23)

其侧向偏向力系数CY有，

$C_{Y} = C_{Y0} + C_{Y}^\beta \beta$ (24)

由上述（23）（24）明显可以得到控制面控制的数学物理解释。

接下来介绍三个重要控制面的作用：副翼（Ailerons），用于控制飞机的滚转运动。通常安装在机翼后缘的外侧部分。升降舵（Elevators）：位于水平尾翼上，用来控制飞机的俯仰运动。方向舵（Rudder）：位于垂直尾翼上，用来控制飞机的偏航运动。

综上，之所以产生了侧向力矩。是因为，控制面在航行方向上始终存在着与机身轴线的恒定夹角，导致产生了偏向力。设这个相对夹角为原角度加撇。

下面，我们将用副翼、升降舵和方向舵三个控制面来分别控制滚转、俯仰或偏航姿态。由式（23）（24）结合力矩的定义$M = F \times r$,得到三个自由度的应控制面偏转而多出来的偏转力矩表达，

$\left\{\begin{matrix}M_{\theta} = \frac{1}{2} \rho A_{a}r_{a}(m_{0}^\theta+ m^\theta \theta' ) \\M_{\beta} = \frac{1}{2} \rho A_{e}r_{e}(m_{0}^\beta+ m^\beta \beta' ) \\M_{\phi} = \frac{1}{2} \rho A_{r}r_{r}(m_{0}^\phi+ m^\phi \phi' ) \end{matrix}\right.$ (25)

下面对上述的全控制系统绘制简单控制框图（图-7），只展示一个自由度，

（图-7 控制系统框图）

图中，r表示预期控制角度、δ表示控制面夹角、e表示控制舵机的PWM信号、α1表示控制面转动角、α2表示该自由度实际偏转角。图中的“动力学模型”即为（25）式。将力矩代入Newton动力学方程，对其进行Laplace变换，假定输入信号为r，先只取俯仰角方向方程，假定输入信号为预期飞行角度r，其输出信号为控制面夹角,并令$K_{\beta} = \frac{1}{2} \rho A_{e} r_{e} m^\beta$,有

$\delta(s_{r}) = \frac{K_{\beta}}{Is_{r}^2}$ (26)

式中sr为经过拉氏变换后的r信号，其中的控制器可以采用基础PID控制器，其传递函数为

$e(\delta) = K_{p} + \frac{K_{i}}{\delta} + K_{d}\delta$ (27)

对于控制舵机可以将其看作一个一阶惯性环节，即

$\alpha_{1} = \frac{1}{\tau e + 1}$ (28)

其中$\tau$为时间系数，其中气动传函为（12）式变形得到，它用来描述实际运动情况。具体操作流程较为复杂，本文在此不做过多赘述。由控制流程化简原理，将（26）（27）（28）式相乘可得，

$\alpha_{1}(s_{r}) = \frac{1}{\tau K_{p} + \frac{\tau K_{i}}{\delta}+\frac{\tau K_{d}K_{\beta }}{I S_{r}}\delta + 1 }$ (29)

若考虑在IMU传感器下形成的反馈闭环回路，设其增益系数为$K_{f}$,其中$\delta_{0}$为初态值，满足线性关系，

$\delta_{r} = K_{f} \alpha_{1} + \delta_{0}$ (30)

作用于方程（29）式，得，

$\alpha_{1}(s_{r}) = \frac{1}{\tau K_{p} + \frac{\tau K_{i}}{\delta}+\frac{\tau K_{d}K_{\beta }}{I (K_{f} \alpha_{1} + \delta_{0})}\delta + 1 }$ (31)

根据（31）式设计迭代算法，就可以初步实现对飞行姿态的控制。

上述的方法可以推广到多个自由度的情况，经过反复的调试和调整控制参数。理论上可以实现初步的控制。尽管实际情况的复杂程度可能超过我们的认知，但在理论设计阶段应侧重于原理和方案的可行性，即“抓主要矛盾，忽略次要因素”。

四 声速\超声速飞行

某些高级爱好者的飞机，有时达到声速甚至更高。那么，对于超声速运动的物体，在和空气接触时可能会产生激波，这部分波的播阻会影响到完整的飞行过程。为了简化模型，我们可以做出以下基本假设：

1.全过程考虑空气为理想气体。满足$PV = nRT$。

2.视流体为可压缩流，即$\frac{\partial \rho }{\partial t} = 0$。

3.假设大气条件符合国际标准大气（ISA）模型。

4.假设流动是定常的，即所有除超音速流体具有的特殊性质参数，其余的流动参数不随时间变化。

考虑飞行器在水平方向以接近1Ma飞行，此时来流总平行于飞行器航向。若此时继续加速，机体头锥处高压区空气将被继续压缩，直到突破音障（Sound barrier）。并激发正激波（normal shock），即波面垂直于来流方向的激波。（如图-8）

（图-8 钝头体的激波产生示意图）

从本质上来讲，激波也存在Doppler效应。此时波源速度明显大于当地声速，即$v_{s} > v_{0}$,

如图-9所示，波源从A点出发沿正方向高速运动。从某一时刻t0发出一道横波。取相同时间间隔T。不难发现，波面间满足后者落后前者$2\pi$个相位差。

（图-9 正激波示意图）

在t时刻时波源位于s处，作各波面的包络直线，可以得到一圆锥面，设α为圆锥半顶角,分析明显可得，

$\sin \alpha = \frac{AP}{AS}=\frac{v}{v_{s}}$ (32)

那么，对于正激波由前面提到的“基本假设”我们可以利用一维定常等熵流动的守恒定律来解决该问题。

由一般可压缩流体质量守恒（连续性）方程，

$\frac{\partial \rho }{\partial t} + \nabla\cdot (\rho v) = 0$ (33)

再考虑N-S（动量守恒）方程，

$\rho \frac{\partial v}{\partial t} + \rho (v \cdot \nabla )v = \rho F -\nabla p + \eta \nabla ^2 v$ （34）

以及能量守恒方程，

$\frac{\partial E}{\partial t} = \nabla \cdot [(E + p) v ] = \nabla \cdot (\tau \cdot v) - \nabla \cdot q + \rho f \cdot v$ (35)

该方程中τ为表面应力张量，其中总能量E的定义为，

$E = \rho e + \frac{1}{2} \rho v^2$ （36）

（36）式中，q为热通量向量（热传导项），e为单位质量内能。对于等熵流动过程，单位质量内能只和温度有关，

$e = \sum h = c_{v}T$

由Mayer关系$c_{p} - c_{v} = R$和比热比的定义，R为气体常数，化简得

$e = \frac{1}{\gamma -1}RT$ (37)

联立（33）（34）（35）（36）（37）式和理想气体方程，得到本构方程组。由实际情况下的边界条件，在考虑动参考系下和飞行器的相对运动的情况下，可使用该方程组对超音速流场中的飞行器进行分析。

在大气中高速运动的飞行器的力学分析和姿态控制的难度巨大，肯定无法用简单的一篇文章就可以一概而论的。本文对此只作基础性说明，并不深入研究。

五 结束语

飞行力学是以空气动力学为基础的研究运动规律和特性的飞行器设计领域的系统、综合性学科，在失速机动战斗机、无人作战飞机、大型运输机、大型客机等领域有着广泛的应用。对于爱好者而言它可以方便我们理解自己制作模型飞机的稳定性和飞行过程。分别我们进行总体设计。

而姿态控制方法则关乎到飞行器的可控性及稳定性，二者密不可分。本文以深受广大爱好者们欢迎的固定翼航模飞行器为切入点，详细研究了其飞行的力学过程，并建立数学模型加以分析。同时指出了控制面控制方案的合理性，最后提及了些许超声速飞行的原理。文章兼具教学作用，旨在为部分初入江湖的爱好者们阐明飞行原理。

航空航天类爱好早期由国外引进，我国的起步相较于某些发达国家较晚。但近年来，随着国民科教素养的提升和政策的放开。我们的航模类活动也在如火如荼的发展着（尤其是年轻一代）。

（图-10 国外某爱好者和他DIY模型的合影）

（图-11 关于低空经济的政府发文）

而且，最近我国也出台了“低空经济”政策，其中“固定翼飞行器”身在其列（文献[3]）。所以，我国爱好者的探索和科创之路依然任重而道远。不过笔者相信在不久的将来，这些曾在我们手里的“玩物”将终有一天会飞向未来，飞向世界。

参考文献

[1] 三水合番.动态滑翔—— 245m/s 的无动力滑翔机[O].2023.https://www.kechuang.org/t/88681

[2] 韩艳铧,陆宇平.变质心控制高速滑翔弹动力学与控制律研究[J].弹道学报, 2012, 24(2):XXXXI:10.3969/XXXXsn.1004-499X.2012.02.003.

[3] 人民日报.首次写入政府工作报告——“低空经济”加速起飞[N].2024.XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

[4] 陈俊谷,刘俊辉,单家元,等.再入滑翔飞行器速度约束减速制导方法[J].兵工学报:0-null[2025-07-05].DOI:10.12382/bgxb.2025.0202.

[5] 李擎,苏中,范军芳.固定推力和总冲约束的高空飞行器姿态控制[J].北京信息科技大学学报：自然科学版, 2012, 27(5):XXXXI:CNKI:SUN:BJGY.0.2012-05-006.

[6]张金学,掌明,李媛媛.基于四元数法的固定翼微型飞行器姿态控制[J].计算机测量与控制, 2012, 20(7):XXXXI:CNKI:SUN:JZCK.0.2012-07-040.

[7] 强元棨.经典力学（上下册）[M].北京:科学出版社.2003,20