数学建模-寻找失事客机黑匣子(2)

问题描述

原问题请参考第一篇博客

假设黑匣子落水之后，不考虑洋流流动对黑匣子沉降过程的影响，建立模型描述黑匣子在水中沉降过程轨迹。如图1所示，假设黑匣子落水点所对应的海底位置为1，落水时沿着图1中指定的虚线方向沉海，给出黑匣子沉在海底的位置，并指出在图形中的哪个区域范围。

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

符号说明

模型建立与求解

问题二

模型分析

现今为了解决水上事故后定位打捞黑匣子困难以及陆地事故后黑匣子存活率达不到100%的难题，抛放式黑匣子应运而生。此类黑匣子可通过其坠毁感知传感器监控飞机事故时触地或坠海瞬间的特征参数异常变化，迅速控制其与机体抛放分离，原理与汽车在撞击瞬间释放安全气囊相似。黑匣子在事故瞬间离机后，如果落在陆地，可避免机体残骸的冲击和火烧等破坏。

所以我们假设黑匣子在飞机坠海瞬间即与飞机脱离，保持飞机坠海前瞬间的速度，水平速度约为 $171m/s$ ,竖直速度约为 $156.35m/s$ 且黑匣子与海面接触后瞬间速度不变。经资料查询假设,黑匣子质量 $m’=50kg$ ,尺寸为 $50cm \times 20cm\times 15cm$ 。即体积为 $V' = 0.015m^3$

黑匣子在海中沉降过程中，受到如下力：

水的粘滞阻力 $f'$ $f^{'}$ ，其大小与速度的平方成正比，方向与前进速度相反，可沿竖直与水平速度方向 $v'_1,v'_2$ $v_{1}^{'}, v_{2}^{'}$ 上分解为 $f_1',f_2'$ $f_{1}^{'}, f_{2}^{'}$
- 水平方向粘滞阻力公式 $f_2'=k'_2 v_2^2$
- 竖直方向粘滞阻力公式 $f_1'=k'_1 v_1^2$
受到水的浮力 $F'$ $F^{'}$ ，其大小恒定，方向始终竖直向上
- $F'=\rho_水 g V=1000*9.8*0.015=147N$
受到自身的重力 $G'$ $G^{'}$ ，其大小恒定，方向始终竖直向下
- $G'=m'g=490N$

黑匣子在水中的受力分析图

经查阅相关资料得出如下系数，

竖直方向粘滞阻力系数 $k'_1=0.2$

水平方向粘滞阻力系数 $k'_2=0.06$

第i时间段末的速度

$v'_1(i),v'_2(i)$ 分别表示第i时间段末的竖直与水平的速度

第i时间段的粘滞阻力

$f'_1(i),f'_2(i)$ 分别表示第i时间段末的竖直与水平的粘滞阻力

第i时间段的竖直加速度

$a_1'(i)=g-\frac{F'+f'_1(i-1)}{m'}$

第i时间段的水平加速度

$a_2'(i)=\frac{f'_2(i-1)}{m'}$

第i时间段末时刻深度，垂直速度，水平速度，水平位移同之前一样

初值说明

初始深度 $h'(0)=0$

初始竖直速度 $v_1'(0)=156.35m/s$

初始水平速度 $v_2'(0)=171m/s$

代码实现

from functools import lru_cache
k1 = 0.2
k2 = 0.06
g = 9.8
m = 50
F = 147
dt = 0.001

@lru_cache
def f1(i):
    return k1*v1(i)**2

@lru_cache
def f2(i):
    return k2*v2(i)**2

@lru_cache
def a1(i):
    return g - (F + f1(i-1)) / m


@lru_cache
def a2(i):
    return f2(i-1) / m


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


@lru_cache
def v2(i):
    if i == 0:
        return 171
    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 0
    ret = h(i-1) - 0.5 * dt * (v1(i) + v1(i-1))
    if ret < -5000:
        return -5000
    return ret


with open('output.csv', 'w') as f:
    H = h(0)
    i = 1
    while H != -5000:
        H = h(i)
        S = s(i)
        V1 = v1(i)
        V2 = v2(i)
        A1 = a1(i)
        A2 = a2(i)
        tup = (i*dt, H, S,
               V2, V1, (V1**2+V2**2)**0.5,
               A2, A1, (A1**2+(A2+g)**2)**0.5,0)
        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