Balistická křivka¶
příklad výpočtu trajektorie hmotného bodu při šikmém vrhu se započtením odporu vzduchu (balistické křivky) numerickým řešením pohybové rovnice
importujeme knihovny numpy (na numerické výpočty) a matplotlib (na kreslení grafů)
In [3]:
import numpy as np
import matplotlib.pyplot as plt
Odporová síla prostředí $F_o = \frac{1}{2}\rho S C_d v^2$ (Newtonův zákon)
In [5]:
def Fo(v,rho,Cd,S):
return 0.5*rho*Cd*S*v**2
konstanty
In [7]:
g = 9.82 #gravitační zrychlení m/s^2
Cd = 0.5 #součinitel odportu pro kouli v oblasti platnosti Newtonova zákona Fo ~ v^2
rho = 1.29 #hustota vzduch (kg/m^3)
mu = 2e-5 #viskozita vzduchu (Pa s)
d = 0.005 #průměr střely (5 mm)
rho_bullet = 11340 #hustota střely z olova (kg/m^3)
S = np.pi*(d/2)**2 #průřez střely
m = 4/3*np.pi*(d/2)**3*rho_bullet #hmotnost střely
v0 = 150 #počáteční rychlost střely (m/s)
alpha = 45*np.pi/180 #úhel výstřelu (45 stupňů) -> převedeno do rad.
Re = rho*d*v0/mu #Reynoldsovo číslo
In [8]:
print(f"hmotnost střely: {m*1000:4.2f} g")
print(f"Reynoldsovo číslo: {Re}")
hmotnost střely: 0.74 g Reynoldsovo číslo: 48375.0
Reynoldsovo číslo odpovídá turbulentnímu proudění v oblasti platnosti Newtonova zákona, tj. $F_o \propto v^2$.
Zvolíme max. čas $t_{max}$, po který budeme pohyb sledovat, a časový krok $dt$ pro numerické řešení.
In [11]:
tmax = 50 #doba, po kterou šikmý vrh sledujeme (zde 50 s)
dt = 0.001 #časový krok pro numerický výpočet derivace (zde tisícina sekundy)
Vytvoříme pole t obsahující časy,
pole x a y obsahující souřadnice v těchto časech
pole vx a vy obsahující rychlosti v těchto časech
In [13]:
t = np.arange(0,tmax,dt) #pole časů
x = np.empty(np.size(t)) #pole x-ových souřadnic
y = np.empty(np.size(t)) #pole y-ových souřadnic
vx = np.empty(np.size(t)) #pole x-ových složek rychlostí
vy = np.empty(np.size(t)) #pole y-ových složek rychlostí
Můžeme ještě vytvořit pole v obsahující velikost rychlosti v jednotlivých časech, $v=\sqrt{v_x^2+v_y^2}$
a pole s obsahující dráhu, kterou hmotný bod urazil, $s=\int_0^t v(t)\,dt$.
In [15]:
v = np.empty(np.size(t)) #pole velikostí rychlostí
s = np.empty(np.size(t)) #pole uražené dráhy
počáteční podmínky:
- hmotný bod je v čase $t = 0$ v počátku, tj. v bodě $[0,0]$
- v čase $t = 0$ má rychlost o velikosti $v_0$ a vektor rychlosti svírá s osou $x$ úhlel $\alpha$, tj. vektor rychlosti v čase $t = 0\,$ je $\,\vec{v}(t=0) = (v_0\,\text{cos}\alpha,v_0\,\text{cos}\alpha)$.
In [17]:
x[0] = 0 # x-ová souřadnice v čase t = 0
y[0] = 0 # y-ová souřadnice v čase t = 0
v[0] = v0 #velikost rychlosti v čase t = 0
vx[0] = v0*np.cos(alpha) - Fo(v0,rho,Cd,S)/m*np.cos(alpha)*dt/2 # x-ová složka rychlosti v čase t = 0 + dt/2
vy[0] = v0*np.sin(alpha)-g*dt/2 - Fo(v0,rho,Cd,S)/m*np.sin(alpha)*dt/2 # y-ová složka rychlosti v čase t = 0 + dt/2
v[0] = np.sqrt(vx[0]**2+vy[0]**2)
s[0] = 0 #uražená dráha v čase t = 0
Nyní numericky vyřešíme pohybovu rovnici šikmého vrhu $$m\vec{a} = \vec{F}_g + \vec{F}_o,$$ což rozepsáno do složek dává \begin{align} \ddot{x} &= -\frac{1}{2}\frac{\rho C_d v^2}{m}\frac{{\dot{x}}}{v}\\ \ddot{y} &= -g-\frac{1}{2}\frac{\rho C_d v^2}{m}\frac{\dot{y}}{v}\\ \end{align}
In [19]:
i = 0 #index pro cyklus přes jednotlivé časy
imax = np.size(t) #max. hodnota indexu (abychom nepřekročily rozsah pole)
istop = 0 #přepínač kontrolující jestli předmět dopadl na zem, pokud dopadne na zem -> istop=1 a cyklus skončí
while istop == 0 and i < np.size(t)-1: #cyklus běžéí dokud předmět nespadne na zem nebo neprojdu celé pole časů
i+=1
x[i] = x[i-1]+vx[i-1]*dt #update x-ové souřadnice pro čas t -> d + dt
y[i] = y[i-1]+vy[i-1]*dt #update y-ové souřadnice pro čas t -> d + dt
if y[i]<0: #spadlo to na zem
istop = 1
imax = i
v[i] = np.sqrt(vx[i-1]**2+vy[i-1]**2) #velikost rychlosti
vx[i] = vx[i-1]-Fo(v[i],rho,Cd,S)/m*vx[i-1]/v[i]*dt #update x-ové složky rychlosti pro čas t -> d + dt
vy[i] = vy[i-1]-g*dt-Fo(v[i],rho,Cd,S)/m*vy[i-1]/v[i]*dt #update y-ové složky rychlosti pro čas t -> d + dt
s[i] = s[i-1]+v[i]*dt #aktuálně uražená dráha
Nakreslíme trajektorii hmotného bodu
In [21]:
plt.plot(x[0:imax],y[0:imax])
plt.xlabel("x")
plt.ylabel("y")
plt.xlim(0,)
plt.ylim(0,)
plt.title("trajektorie")
plt.show()
Pro srovnání můžeme vypočítat, jak by vypadala trajektorie bez odporu vzduchu
In [23]:
x_b = np.empty(np.size(t)) #pole x-ových souřadnic bez odporu vzduchu
y_b = np.empty(np.size(t)) #pole y-ových souřadnic bez odporu vzduchu
vx_b = np.empty(np.size(t)) #pole x-ových složek rychlostí bez odporu vzduchu
vy_b = np.empty(np.size(t)) #pole y-ových složek rychlostí bez odporu vzduchu
v_b = np.empty(np.size(t)) #pole velikostí rychlostí bez odporu vzduchu
s_b = np.empty(np.size(t)) #pole uražené dráhy bez odporu vzduchu
x_b[0] = 0 # x-ová souřadnice v čase t = 0 bez odporu vzduchu
y_b[0] = 0 # y-ová souřadnice v čase t = 0 bez odporu vzduchu
v_b[0] = v0 #velikost rychlosti v čase t = 0 bez odporu vzduchu
vx_b[0] = v0*np.cos(alpha) # x-ová složka rychlosti v čase t = 0 + dt/2 bez odporu vzduchu
vy_b[0] = v0*np.sin(alpha)-g*dt/2 # y-ová složka rychlosti v čase t = 0 + dt/2 bez odporu vzduchu
v_b[0] = np.sqrt(vx_b[0]**2+vy_b[0]**2) #velikost rychlosti v čase t = 0 bez odporu vzduchu
s_b[0] = 0 #uražená dráha v čase t = 0 bez odporu vzduchu
i = 0 #index pro cyklus přes jednotlivé časy
imax_b = np.size(t) #max. hodnota indexu (abychom nepřekročily rozsah pole)
istop = 0 #přepínač kontrolující jestli předmět dopadl na zem, pokud dopadne na zem -> istop=1 a cyklus skončí
while istop == 0 and i < np.size(t)-1: #cyklus běžéí dokud předmět nespadne na zem nebo neprojdu celé pole časů
i+=1
x_b[i] = x_b[i-1]+vx_b[i-1]*dt #update x-ové souřadnice pro čas t -> d + dt
y_b[i] = y_b[i-1]+vy_b[i-1]*dt #update y-ové souřadnice pro čas t -> d + dt
if y_b[i]<0: #spadlo to na zem
istop = 1
imax_b = i
v_b[i] = np.sqrt(vx_b[i-1]**2+vy_b[i-1]**2) #velikost rychlosti
vx_b[i] = vx_b[i-1] #update x-ové složky rychlosti pro čas t -> d + dt
vy_b[i] = vy_b[i-1]-g*dt #update y-ové složky rychlosti pro čas t -> d + dt
s_b[i] = s_b[i-1]+v_b[i]*dt #aktuálně uražená dráha
Graf srovnání trajektorie bez odporu vzduchu a se započítaným odporem vzduchu
In [25]:
plt.plot(x_b[0:imax_b],y_b[0:imax_b],c='blue',ls='dashed',label="bez odporu vzduchu")
plt.plot(x[0:imax],y[0:imax],c='blue', label="s odporem vzduchu")
plt.xlabel("x")
plt.ylabel("y")
plt.xlim(0,)
plt.ylim(0,)
plt.title("trajektorie")
plt.legend()
plt.show()
print(f"Vzdálenost místa dopadu: {x[imax]:5.1f} m (se započítaným odporem vzduchu)")
print(f"Vzdálenost místa dopadu: {x_b[imax_b]:5.1f} m (bez odporu vzduchu)")
Vzdálenost místa dopadu: 221.7 m (se započítaným odporem vzduchu) Vzdálenost místa dopadu: 2291.3 m (bez odporu vzduchu)
Časová závislost velikosti rychlosti
In [27]:
plt.plot(t[0:imax_b],v_b[0:imax_b],c='green',ls='dashed',label="bez odporu vzduchu")
plt.plot(t[0:imax],v[0:imax],c='green',label="s odporem vzduchu")
plt.xlabel("čas (s)")
plt.ylabel("rychlost (m/s)")
plt.xlim(0,)
plt.ylim(0,)
plt.title("velikost rychlosti")
plt.legend()
plt.show()
print(f"čas dopadu : {t[imax-1]:4.2f} s (se započítaným odporem vzduchu)")
print(f"čas dopadu : {t[imax_b-1]:4.2f} s (bez odporu vzduchu)")
čas dopadu : 9.10 s (se započítaným odporem vzduchu) čas dopadu : 21.60 s (bez odporu vzduchu)
uražená dráha
In [29]:
plt.plot(t[0:imax_b],s_b[0:imax_b],c='red',ls='dashed',label="bez odporu vzduchu")
plt.plot(t[0:imax],s[0:imax],c='red',label="s odporem vzduchu")
plt.xlabel("čas (s)")
plt.ylabel("dráha (m)")
plt.xlim(0,)
plt.ylim(0,)
plt.title("uražená dráha")
plt.legend()
plt.show()
print(f"celková uražená dráha : {s[imax-1]:4.2f} m (se započítaným odporem vzduchu)")
print(f"celková uražená dráha : {s_b[imax-1]:4.2f} m (bez odporu vzduchu)")
celková uražená dráha : 324.26 m (se započítaným odporem vzduchu) celková uražená dráha : 1133.99 m (bez odporu vzduchu)