[카테고리:] 물리

투사체 질량  m
항력  F=-C_d v=-C_d \sqrt{\dot{x}+\dot{y} }
초기속도  v=v_0, θ=θ_0

m \ddot{x}+C_d \dot{x} =0
m \ddot{y}+C_d \dot{y} =-mg

m \ddot{x}+C_d \dot{x} =0

특성방정식 m \lambda ^2+C_d \lambda =0
\lambda =0, -C_d/m

x(t)=c_1 e^0+c_2 e^{- \frac{C_d}{m} t}=c_1+c_2 e^{- \frac{C_d}{m} t}

초기조건
x'(t)=c_2 {- \frac{C_d}{m}} e^{-\frac{C_d}{m} t}
x'(0)=v_0  cos⁡(θ_0), x_0=0

c_2=- \frac{m}{C_d}  v_0  \cos⁡{θ_0}, c_1=\frac{m}{C_d}  v_0  \cos⁡{θ_0} x(t)=\frac{m}{C_d}  v_0  \cos⁡{θ_0 } - \frac{m}{C_d}  v_0  \cos⁡{θ_0 } \cdot e^{- \frac{C_d}{m} t} my C_d \ddot{y}=-mg=r(t)

r(t)=0으로 놓고 제차 상미방으로 풀이

특성방정식 m \lambda^2+C_d \lambda=0

\lambda=0, -\frac{C_d}{m} y_h (t)=c_1 e^0+c_2 e^{-\frac{C_d}{m} t}=c_1+c_2 e^(-\frac{C_d}{m} t) y_h (t)=\frac{m}{C_d}  v_0  \sin{θ_0}-(\frac{m}{C_d}  v_0  \sin⁡{θ_0}+(m^2 g)/C_d ) \cdot e^{-\frac{C_d}{m} t}

r(t)=-mg=k꼴이므로 y_p (t)=A로 놓으면 답이 안 나옴
y_p (t)=Ax+B라고 하자

C_d \cdot (A)=-mg, B=\frac{m^2 g}{{C_d}^2}

y_h (t)=\frac{m}{C_d}  v_0  \sin⁡{θ_0}- (\frac{m}{C_d}  v_0  \sin⁡{θ_0}+\frac{m^2 g}{C_d} ) \cdot e^{-\frac{C_d}{m} t}
y_p (t)=-\frac{mg}{C_d}  t+\frac{m^2 g}{{C_d}^2}

y(t)=\frac{m}{C_d}  v_0  \sin⁡(θ_0)-(\frac{m}{C_d}  v_0  \sin(θ_0)+(m^2 g)/C_d ) \cdot e^{-\frac{C_d}{m} t} -\frac{mg}{C_d}  t+\frac{m^2 g}{{C_d}^2}

import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)

v0 = 100.0          # m/s
#a0                 # rad
m = 1.0             # kg
#Cd                 # kg/s
g = 9.8             # m/s^2

cond = [(1/3 * math.pi, 0.1), (1/4 * math.pi, 0.1), (1/6 * math.pi, 0.1)]   # (a0, Cd)

dt = 0.001          # s

def doit(a0, Cd):
    t = [0.0]
    x = [0.0]
    y = [0.0]

    while y[-1] >= 0.0:
        t.append(t[-1] + dt)
        ex = math.exp(-Cd / m * t[-1])
        x.append(
            m / Cd * v0 * math.cos(a0)
            - m / Cd * v0 * math.cos(a0) * ex
        )
        y.append(
            m / Cd * v0 * math.sin(a0)
                 - (m / Cd * v0 * math.sin(a0) + m**2 * g / Cd**2) * ex
                 - m * g / Cd * t[-1]
                 + m**2 * g / Cd**2
                 )
    return (x, y)

leg = []

fig, ax = plt.subplots()

for a0, Cd in cond:
    (x, y) = doit(a0, Cd)
    ax.plot(x, y)
    leg.append(f"v(0) = {v0:.1f} [m/s], θ(0) = {a0/math.pi*180.0:.1f} [deg], Cd = {Cd:.1e} [kg/s], m = {m:.2e} [kg]")

plt.xlim(0.0, v0**2 / g)
plt.ylim(0.0, v0**2 / g / 2)
plt.xlabel("x [m]")
plt.ylabel("y [m]")
plt.gca().set_aspect('equal', 'box')
plt.legend(leg, fontsize="x-small")
ax.xaxis.set_major_locator(MultipleLocator(100))
ax.xaxis.set_major_formatter('{x:.0f}')
ax.xaxis.set_minor_locator(MultipleLocator(10))
ax.yaxis.set_major_locator(MultipleLocator(100))
ax.yaxis.set_major_formatter('{x:.0f}')
ax.yaxis.set_minor_locator(MultipleLocator(10))
plt.savefig('posplot.svg', transparent=True)

import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)

v0 = 100.0          # m/s
a0 = 5/12 * math.pi  # rad
m = 0.15             # kg
Cd = 0.1            # kg/s
g = 9.8             # m/s^2

dt = 0.001          # s

def doit(a0, Cd):
    t = [0.0]
    vx = [0.0]
    vy = [0.0]
    y = 0.0

    while y >= 0.0:
        t.append(t[-1] + dt)
        vx.append(
            v0 * math.cos(a0) * math.exp(-Cd / m * t[-1])
        )
        vy.append(abs(
                 (v0 * math.sin(a0) + m * g / Cd) * math.exp(-Cd / m * t[-1])
                 - m * g / Cd
                 ))
        y = m / Cd * v0 * math.sin(a0) \
                 - (m / Cd * v0 * math.sin(a0) + m**2 * g / Cd**2) * math.exp(-Cd / m * t[-1]) \
                 - m * g / Cd * t[-1] \
                 + m**2 * g / Cd**2
    return (t, vx, vy)

leg = []

fig, ax = plt.subplots()

(t, vx, vy) = doit(a0, Cd)
ax.plot(t, vx)
ax.plot(t, vy)
vt = []
for va, vb in zip(vx, vy):
    vt.append(math.sqrt(va ** 2 + vb ** 2))
ax.plot(t, vt)
leg.append(f"vx, v(0) = {v0:.1f} [m/s], θ(0) = {a0/math.pi*180.0:.1f} [deg], Cd = {Cd:.1e} [kg/s], m = {m:.2e} [kg]")
leg.append(f"vy, v(0) = {v0:.1f} [m/s], θ(0) = {a0/math.pi*180.0:.1f} [deg], Cd = {Cd:.1e} [kg/s], m = {m:.2e} [kg]")
leg.append(f"v, v(0) = {v0:.1f} [m/s], θ(0) = {a0/math.pi*180.0:.1f} [deg], Cd = {Cd:.1e} [kg/s], m = {m:.2e} [kg]")

plt.xlim(0.0, t[-1])
plt.ylim(0.0, v0)
plt.xlabel("t [s]")
plt.ylabel("v [m/s]")
#plt.gca().set_aspect('equal', 'box')
plt.legend(leg, fontsize="x-small")
ax.xaxis.set_major_locator(MultipleLocator(1))
ax.xaxis.set_major_formatter('{x:.0f}')
ax.xaxis.set_minor_locator(MultipleLocator(0.1))
ax.yaxis.set_major_locator(MultipleLocator(10))
ax.yaxis.set_major_formatter('{x:.0f}')
ax.yaxis.set_minor_locator(MultipleLocator(1))
plt.savefig('velplot2.svg', transparent=True)

\text{Constants}\\ C_d=0.47\, \text{공의 항력계수}\\ \rho=1.293 [kgm^{-3}]\, \text{공기의 밀도, 1기압, 섭씨 25도}\\ D=8730 [kgm^{-3}]\, \text{Brass, 섭씨 0도}// \text{Physical Variables}\\ L_t=0.7 [m]\, \text{잡아당기는 거리}\\ t_r=0.15 [s]\, \text{반응속도, 프로게이머 기준}\\ V_c=100 [ms^{-1}]\, \text{눈을 관통하기 위한 접촉 시점의 요구 속력}\\ r_e=0.012 [m]\, \text{안구의 반지름}\\ r_b=r_e\, \text{공의 반지름}

\text{Equations}\\ \text{새총 방정식}\\ \varepsilon =\frac{L}{2} \; \therefore \frac{d\varepsilon}{\varepsilon}=\frac{x}{L}=\frac{dx}{L} \quad d\varepsilon=\frac{\varepsilon}{L} dx \qquad \dot{\varepsilon} = \frac{\varepsilon}{L} \dot{x} \\ E^k_b=\frac{1}{2} M_b \dot{x}^2 \\ E^k_s=\frac{1}{2} \int_0^L \sigma d\varepsilon \left( \frac{\varepsilon}{L} \dot{x} \right) ^2 \, dx=\frac{1}{6} \sigma L \dot{x}^2=\frac{1}{6}M_b \dot{x}^2\\ E^k = \frac{1}{2} \left( M_b + \frac{1}{3} M_s \right) \dot{x}^2\\ E^p=\frac{1}{2}kx^2\\ E=\frac{1}{2} \left( M_b + \frac{1}{3} M_s \right) \dot{x}^2 + \frac{1}{2} k x^2\\ \frac{dE}{dx} = 0 =\frac{d}{dx} \left\{ \left( \frac{1}{2} M_b+\frac{1}{6} M_s \right) \dot{x}^2 + \frac{1}{2} k x^2 \right\} = \left( M_b+ \frac{1}{3} M_s \right) \ddot{x} + kx \qquad \ddot{x}=\frac{kx}{M_b + \frac{1}{3} M_s } \\ v = \int_0^t \ddot{x} \, dx = \frac{k}{M_b+\frac{1}{3} M_s} \cdot \frac{1}{2} t^2 \\ \text{항력 방정식}\\ F_A=-\frac{1}{2}C_d \rho A_b v^2\\ a_A= -\frac{-3 C_d \rho }{8Dr} v^2\\ \dot{x} = \int_0^T - \frac{3 C_d \rho }{8Dr} \dot{x}^2 t dt + v_0 =-\frac{3 C_d \rho }{8Dr} \dot{x}^2 +v_0 \\ v_0 = \frac{3 C_d \rho }{8Dr} v^2 + v