Re: 제로부터 시작하는 2차원 새총역학

투사체 질량  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)