阻力模拟？

发表于 更新于 分类于 本文字数： 2.1k 阅读时长 ≈ 2 分钟

物理作业

比较炮弹的实际弹道与理想抛物线的差异

推导

假设质点只受重力、阻力

重力加速度取\(|\vec{g}|=9.8m/{s^2}\) 质点质量为\(m\quad\) 空气阻力\(f\)与速度\(v\)的关系

\[f=\frac{1}{2}C\rho S{ {v}^2}\]

令\(v^2\)的系数\(\frac{1}{2}C\rho S=k\)

有

\[ \vec{G}=m\vec{g} \]

\[ \scriptsize \begin{cases}|\vec{f}|&=&k{v^2}\\\frac{\vec{f} }{|\vec{f}|}&=&-\frac{\vec{v}}{|\vec{v}|}\end{cases} \]

则对于时刻\(t_i\)，令其经历极小时间段\(\Delta t\)后得到时刻为\(t_{i+1}\)，计算时刻\(t_{i+1}\)有

\[ \vec{G_i}=\vec{G_{i+1} }=m\vec{g} \]

\[ \scriptsize \begin{cases} |\vec{f_{i+1} }|&=&k{v_{i}^2}\\ \frac{\vec{f_{i+1} } }{|\vec{f_{i+1} }|}&=&-\frac{\vec{v_i} }{|\vec{v_i}|} \end{cases} \]

\[ \vec{a_{i+1} }=\frac{\vec{ {F_{合} }_ {i+1} } }{m}=\frac{\vec{G_{i+1} }+\vec{f_{i+1} } }{m} \]

\[ \vec{v_{i+1} }=\vec{v_{i} }+\vec{a_{i+1} }\cdot \Delta t \]

\[ \vec{x_{i+1} }=\vec{x_i}+\vec{v_{i+1} }\cdot \Delta t \]

令\(t_0\)时 \[ \scriptsize \begin{cases} t_0&=&0\\ \vec{x_0}&=& \begin{bmatrix} 0\\ 0 \end{bmatrix}\\ \vec{v_0}&=& \begin{bmatrix} v_x\\ v_y \end{bmatrix}\\ \vec{ {F_{合} }_0}&=&\vec{G}+\vec{f}=m\vec{g}\\ \vec{a_0}&=& \begin{bmatrix} 0\\ |\vec{g}| \end{bmatrix} \end{cases} \]

实践

Python程序

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
import matplotlib.pyplot as plt
import numpy as np

t=1

print("v=(x,y)=",end='')
v=input()[1:-1].split(',')
v=np.array([float(v[0]),float(v[1])])

print("\nm=",end='')
m=float(input())

print("\nk=",end='')
k=float(input())

print("\nstep=",end='')
step=int(input())

def draw (flag):

    vv=v

    x=np.array([0,0],ndmin=2)

    g=9.8

    a=np.array([[0,-g],[0,0]])

    for i in range(1,step+1):
        if flag:
            a[1]=-vv*(k**2)/m
        vv=vv+a.sum(axis=0)*t
        x=np.concatenate((x,[x[i-1]+vv*t]),axis=0)

    x=np.array(x)

    plt.plot(x[...,0],x[...,1])

draw(1)
draw(0)

plt.show()

输入

1
2
3
4
(1,2)
1
0.1
100
1
2
3
4
5
6
7
v=(x,y)=(1,2)

m=1

k=0.1

step=100

\[ \scriptsize \begin{cases} \vec{v_0}= \begin{bmatrix} v_x\\ v_y \end{bmatrix}= \begin{bmatrix} 1m/s\\ 2m/s \end{bmatrix} \\ m=1kg\\ k=0.1kg/m\\ 缺省\Delta t=1s \end{cases} \]

计算至\(i\le step\)，即

\[ i\in \{i|0\le i\le step,i\in N\} \]

时所有的\(\vec{x_i}\)并绘图如下：

missile

其中蓝色曲线为有阻力情况，橙色曲线为无阻力情况

可以明显看出，等时间时，有阻力飞行距离更近、下落更快