问题描述

假设有一架飞机在高空中飞行时突然发生事故，此时飞行高度为10000米，飞行速度是800公里/小时，航向东北方向45°，飞机在地面的投影位置为南纬22.0度，东经88.0度。
请建立模型求解以下问题：

  1. 假定飞机在发生事故时突然失去动力，考虑飞机在降落过程中受到空气气流的影响，建立数学模型，描述飞机坠落轨迹并推测黑匣子的落水点。
  2. 假设黑匣子落水之后，不考虑洋流流动对黑匣子沉降过程的影响，建立模型描述黑匣子在水中沉降过程轨迹。如图1所示，假设黑匣子落水点所对应的海底位置为1，落水时沿着图1中指定的虚线方向沉海，给出黑匣子沉在海底的位置，并指出在图形中的哪个区域范围。

图1 黑匣子在水下沉降过程中的海底剖面图

  1. 考虑洋流流动对黑匣子在水中沉降的影响，建立模型描述在有洋流流动的情况下黑匣子沉降轨迹，并求解出黑匣子沉入水下1000m，2000m和3000m时离落水点的方位。

模型假设

  1. 假设飞机在坠降过程中受飞行员控制机身总保持水平，且整个过程中未受到鸟类撞击等影响，飞机没有并发故障；
  2. 假设飞机坠降中机身温度不发生改变；
  3. 假设飞机在落到水面前保持机身完整，且接触水面瞬间解体，忽略解体时其他物体撞击或爆炸对黑匣子影响；
  4. 假设黑匣子在下落过程中不发生翻滚；
  5. 假设黑匣子在水中沉降不会受到海洋生物影响；
  6. 忽略东南季风对飞机的影响；
  7. 假设客机为波音747-400型客机[3]，且正常载员载货；
  8. 认为飞机机师经验丰富，遇到此类突发事故，会沉着冷静应对，且忽略机上人员、货物因突发事故引起的举动对飞机的影响。

符号说明

  • 重力加速度g=9.8m/s2g=9.8m/s^2

  • 空气摩擦阻力ff

    • 飞机竖直方向向上的摩擦力为f1f_1
    • 飞机水平方向向后的摩擦力为f2f_2

  • 飞机的升力FF

  • 空气密度ρ=1.29kg/m3\rho=1.29kg/m^3

  • 飞机质量m=362875kgm=362875kg

  • 飞机的重力G=mg=3556175NG=mg=3556175N

  • 升力系数CiC_i

  • 竖直方向的飞机速度v1v_1

    • 约定使用v1(t)v_1(t)表示t时刻的飞机竖直速度
    • 已知v1(0)=0v_1(0)=0

  • 水平方向的飞机速度v2v_2

    • 约定使用v2(t)v_2(t)表示t时刻的飞机竖直速度
    • 已知v2(0)=800km/h=222.2m/sv_2(0)=800km/h=222.2m/s

  • 飞机机翼面积S1=541.2m2S_1=541.2m^2

  • 飞机正面迎风面积S2=75m2S_2=75m^2

  • 飞机底面积S3=1000m2S_3=1000m^2

模型建立与求解

问题一

飞机在高空匀速水平飞行时，突发故障，在坠降过程中，飞机速度为v，飞机受到重力G，升力F，摩擦阻力f。假设飞机的到落水的加速度全程未达到最终的解体加速度极限。

  • 飞机的速度可分解为水平方向前进的速度v2v_2，和竖直方向降落的速度v1v_1
  • 重力G方向竖直向下，大小恒定
  • 升力F方向竖直向上，大小变化，见下面的详细分析
  • 摩擦阻力ff可分解成水平和竖直的两个摩擦力。根据水平前进速度v2v_2得到方向向后的水平摩擦力f2f_2；根据竖直的降落速度v1v_1得到方向向上的竖直摩擦力f1f_1

升力计算公式

L=12ρv2SCiL=\frac12 \rho v^2 S C_i

  • ρ\rho为空气密度，标准取值为ρ=1.29kg/m3\rho = 1.29kg/m^3
  • vv为水平方向的速度
  • SS为机翼面积
  • CiC_i为升力系数
计算竖直方向升力L

​ 初始状态下，飞机水平匀速飞行，故竖直方向上无速度，即也无竖直的摩擦阻力f1f_1。竖直方向无速度，故升力与重力二力平衡L=GL=G，即有如下等式成立

12ρv2SCi=mg\frac12 \rho v^2 S C_i=mg

​ 已知空气密度ρ=1.29kg/m3\rho = 1.29kg/m^3，初始的水平速度v2(0)=222.2m/sv_2(0)=222.2m/s，机翼面积S=S1=541.2m2S=S_1=541.2m^2，飞机质量m=362875kgm=362875kg，重力加速度g=9.8m/s2g=9.8m/s^2，未知数仅剩CiC_i，即Ci=2mgρv2S=0.2C_i=\frac{2mg}{\rho v^2 S}=0.2

​ 此时仅剩水平速度v2v_2随着时间而改变，其他量均可完成计算，故该飞机的升力最终计算公式为L=12∗1.29∗541.2∗0.2∗v2=69.8v2L=\frac12 * 1.29 * 541.2 * 0.2 * v^2=69.8v^2，记系数k1=69.8k_1=69.8，则升力计算公式为L=k1v22L=k_1 v_2^2

空气阻力计算公式

F=12ρv2SCdF=\frac12 \rho v^2 S C_d

其中，

  • CdC_d为空气阻力系数，这里取0.08
  • ρ\rho为空气密度
  • SS物体迎风面积
  • vv为物体与空气的相对运动速度。
计算竖直方向空气阻力f1

​ 竖直方向的迎风面积为底面积S3=1000m2S_3=1000m^2

​ 竖直方向的速度为v1v_1

​ 仿照计算升力的步骤，记k2=12∗1.29∗1000∗0.08=51.6k_2=\frac12 * 1.29 * 1000 * 0.08=51.6

​ 则f1=k2v12f_1=k_2 v_1^2

计算水平方向空气阻力f2

​ 水平方向的迎风面积为正面面积S2=75m2S_2=75m^2

​ 水平方向的速度为v2v_2

​ 同理可记，k3=12∗1.29∗75∗0.08=3.87k_3=\frac12 * 1.29 * 75 * 0.08=3.87

​ 则f2=k3v22f_2=k_3v_2^2

坠降过程中速度分解图与受力分析图

微元分析

​ 采用微元法，将整个过程分成n段以Δt=0.001s\Delta t=0.001s为间隔的小段，将每一小段的水平和竖直方向各看成匀加速直线运动，现对其中任意的第i小段。

第i时间段内竖直方向的加速度a1(i)

竖直方向在坠江过程受到向下的重力G，向上的升力F，向上的空气阻力f1f_1，即

在竖直方向上坠降时具有向下的加速度a1a_1，做加速运动即，合力向下，则有

F合=G−f1−LF_合=G-f_1-L

ma1=mg−k1v22−k2v12ma_1=mg-k_1 v_2^2-k_2v_1^2

a1=g−(k1v22+k2v12)ma_1=g-\frac{(k_1 v_2^2 + k_2 v_1^2)}{m}

第i时间段内水平方向的加速度a2(i)

水平方向仅受到向后的空气阻力f2f_2，故其加速度向后，与前进速度相反，做减速运动

F合=f2F_合=f_2

ma2=k3v22ma_2=k_3 v_2^2

a2=k3v22ma_2 = \frac{k_3 v_2^2}{m}

第i时间段末时刻的高度

下一时刻的高度应当是上一时刻的高度减去这一时间段降落的高度

  • h(i)表示第i个时间段末时刻的高度

  • v1(i)v_1(i)表示第i段时间末的物理量v1v_1的值

  • v‾=v1(i)+v1(i−1)2\overline v = \frac{v_1(i)+v_1(i-1)}{2}表示第i段时间内的平均速度

  • 第i段时间内下降的高度Δh=v‾Δt\Delta h=\overline v \Delta t

即第i时间段的末时刻高度为

h(i)=h(i−1)−v1(i)+v1(i−1)2Δth(i)=h(i-1)-\frac{v_1(i)+v_1(i-1)}{2} \Delta t

第i时间段末时刻竖直方向的速度

v1(i)=v1(i−1)+a1(i)∗Δtv_1(i)=v_1(i-1)+a_1(i)*\Delta t

第i时间段末时刻水平方向的速度

v2(i)=v2(i−1)+a2(i)∗Δtv_2(i)=v_2(i-1)+a_2(i)*\Delta t

第i时间段末时刻的水平位移

s(i)=s(i−1)+v2(i)+v2(i)2Δts(i)=s(i-1)+\frac{v_2(i)+v_2(i)}{2} \Delta t

初值说明
  • 质量m=362875kgm=362875kg
  • 重力加速度g=9.8m/s2g=9.8m/s^2
  • 初始高度h(0)=10000mh(0)=10000m
  • 初始竖直方向速度v1(0)=0v_1(0)=0
  • 初始水平方向速度v2(0)=222.2m/sv_2(0)=222.2m/s
  • 初始水平位移s(0)=0s(0)=0

代码编写

初步代码

根据微元法分析得到的几个迭代公式，编写递归代码如下，

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
k1 = 69.8
k2 = 51.6
k3 = 3.87
g = 9.8
m = 362875
dt = 0.001

def a1(i):
    return g - (k1*v2(i)**2 + k2*v1(i-1)**2) / m

def a2(i):
    return (k3*v2(i-1)**2) / m

def v1(i):
    if i == 0:
        return 0
    return v1(i-1) + a1(i) * dt

def v2(i):
    if i == 0:
        return 222.2
    return v2(i-1) - a2(i) * dt

def s(i):
    if i == 0:
        return 0
    return s(i-1) + 0.5 * dt * (v2(i) + v2(i-1))

def h(i):
    if i == 0:
        return 10000
    ret = h(i-1) - 0.5 * dt * (v1(i) + v1(i-1))
    if ret < 0:
        return 0
    return ret

若想要计算任何一个物理量在某个时间段的值，可直接调用上述函数，

这里的递归代码写起来很简单，直接照搬建模结果，但是思考起来较为复杂，如下列计算h(1)的步骤，

假设要求h(1)，那么就要求h(0)，v1(0)，v1(1)

要求v1(1)，那么就要求v1(0)，a1(1)

要求a1(1)，那么就要求v2(1)，v1(0)

要求v2(1)，那么就要求v2(0)，a2(1)

要求a2(1)，那么就要求v2(0)，

此时v2(0)已知，a2(1)便可求解

v2(1)便可求解，a1(1)便可求解

v1(1)便可求解，h(1)便可求解

实际代码中，若不加以优化直接套用上述的递归计算，很容易便出现超过递归最大深度的限制。

并且还会出现重复计算，导致时间复杂度相当高，几乎无法在合适的时间内得到结果。

算法优化

可以使用记忆化搜索的方式来实现对某次计算结果的缓存，来实现提高计算效率。在Python中，使用lru_cache缓存装饰器即可实现函数计算的缓存。

优化后的代码如下，效率大大提升，

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
from functools import lru_cache
k1 = 69.8
k2 = 51.6
k3 = 3.87
g = 9.8
m = 362875
dt = 0.001


@lru_cache
def L(i):
    return k1*v2(i)**2


@lru_cache
def a1(i):
    return g - (L(i) + k2*v1(i-1)**2) / m


@lru_cache
def a2(i):
    return (k3 * v2(i-1) ** 2) / m


@lru_cache
def v1(i):
    if i == 0:
        return 0
    return v1(i-1) + a1(i) * dt


@lru_cache
def v2(i):
    if i == 0:
        return 222.2
    return v2(i-1) - a2(i) * dt


@lru_cache
def s(i):
    if i == 0:
        return 0
    return s(i-1) + 0.5 * dt * (v2(i) + v2(i-1))


@lru_cache
def h(i):
    if i == 0:
        return 10000
    ret = h(i-1) - 0.5 * dt * (v1(i) + v1(i-1))
    if ret < 0:
        return 0
    return ret


with open('output.csv', 'w') as f:
    H = h(0)
    i = 1
    while H != 0:
        H = h(i)
        S = s(i)
        V1 = v1(i)
        V2 = v2(i)
        A1 = a1(i)
        A2 = a2(i)
        # A1 = 0
        # A2 = 0
        tup = (i*dt, H, S,
               V2, V1, (V1**2+V2**2)**0.5,
               A2, A1, (A1**2+(A2+g)**2)**0.5,
               L(i)/(m*g))
        f.write("%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n" % tup)
        print("时间: %fs, 高度: %fm, 水平位移: %fm，水平速度: %f，竖直速度: %f, 合速度：%f，"
        "水平加速度: %f, 竖直加速度: %f，合加速度: %f，G:%f" % tup)
        i += 1

如下图，代码开始进行迭代计算，

image-20220112190505381

当下降至高度为0时，停止计算

image-20220112190946760

数据均保存在output.csv文件中

计算结果
  • 飞机在126.6s时接触水面
  • 整个迫降过程的水平位移为24601.40m
  • 接触水面时，水平速度为170.92m/s
  • 接触水面时，竖直速度为156.35m/s
  • 接触水面时，合速度为231.64m/s

数据分析

轨迹拟合

在excal中，绘制出水平位移与垂直高度的散点图(由于有十几万个点，故相当密集，几乎是一条线)

image-20220112191421209

使用趋势线拟合，如下

image-20220112191727009

使用二次多项式拟合，相关系数达到0.9997，已经可以很好的代表其坠落轨迹

image-20220112191823409

使用四次多项式拟合，相关系数更是直接达到1

地点预测

飞机在坠落之前航向为东北方向45°45\degree，投影位置在南纬22.0°22.0\degree，东经88.0°88.0\degree，理想情况下，飞机在迫降过程的水平方向上仍保持原来的方向。经计算得到飞机的水平位移为24601.40m

飞机将在经线和纬线方向上飞行24601.40/2=17398.44m24601.40/\sqrt{2}=17398.44m

在经线上，纬度每差1度,实地距离大约为111千米

故飞机大约偏离原纬度17.393/111≈0.157°17.393/111\approx 0.157\degree，坠落目标纬度为南纬22−0.157=21.84°22-0.157=21.84\degree

在纬线上，经度每差1度,实际距离为111×cosθ千米，其中θ表示该纬线的纬度，该纬线为南纬22度，故实际距离大约102.92km

故飞机大约偏离原经度17.393/102.92≈0.169°17.393/102.92\approx 0.169\degree，坠落目标经度为东经88+0.169=88.17°88+0.169=88.17\degree

坠落地点如下图所示，

坠落地点

解体判定

​ 如撞击，爆炸等所造成的空中解体是指由于外力造成的。而超过最大过载导致空中解体则是指飞行器在进行机动飞行（非直线运动）时所产生的离心加速度与自身重力加速度的矢量和，单位是“G”。当飞行器在高度速度方向都不变的情况下飞行时，过载是1G。当它开始做变速运动时，升力大于重力的倍数就是多少G。假设一架飞机在做机动动作时，过载达到4G时，这就意味着飞机要承受其重量4倍的力。

​ 飞行器在进行剧烈的机动动作，如翻筋斗、倒飞、高速小半径转弯、极限加速等。所产生的过载超过了机体所能承受的最大过载（客机通常最大是3.5G，战斗机通常最大是9G），就会因为机体结构无法承受如此之大的过载，导致空中解体。

​ 根据计算数据，绘制出客机在坠落过程中升力与重力的比值与时间的曲线如下图所示，

​ image-20220112214938156

​ 如图所示，始终小于3.5G故全程无解体，与假设相一致。

问题二

在另一篇博客

问题三

在另一篇博客

参考文献

空中解体

升力

空气阻力