외부 탄도학 모델링: 해석 및 Python

외부 포탄 운동의 해석적 및 수치 모델링: 세 가지 핵심 사례

외부 포탄 운동은 중력 하에서 수평선과 각도를 이루며 초기 속도로 발사된 물체의 평면 운동을 설명합니다. 세 가지 경우를 고려합니다: 공기 저항 없음, 스토크스 법칙에 따른 선형 저항, 뉴턴 법칙에 따른 제곱형 저항. 해석적 해법은 현실적인 시뮬레이션을 위해 수치 방법과 함께 사용됩니다.

첫 번째 경우—저항 없음—궤도는 포물선입니다. v₀ = 20 m/s, α = 30°일 때, 사거리와 최고 고도는 간단한 공식으로 결정됩니다:

x(t) = v₀ cos(α) t

y(t) = v₀ sin(α) t - (g t²)/2

선형 저항: 스토크스 법칙

선형 저항이 있는 경우, 힘은 속도에 비례합니다: F = -k v. 미분 방정식은 다음과 같습니다:

m dv_x/dt = -k v_x

m dv_y/dt = -m g - k v_y

해결 결과는 지수 함수 형태의 속도 의존성을 보입니다:

v_x(t) = v₀ cos(α) e^(-kt/m)

v_y(t) = (v₀ sin(α) + (m g)/k) e^(-kt/m) - (m g)/k

이들을 적분하면 y(x)를 얻을 수 있습니다. v₀ = 70 m/s, α = 55°, k = 0.04 kg/s, m = 0.5 kg일 경우, 이상적인 경우보다 사거리가 줄어든 비대칭 궤도가 나타납니다.

제곱형 저항: 뉴턴 법칙

고속에서는 제곱형 저항이 주도적입니다: F = -k v². 상미분방정식계는 다음과 같습니다:

dv_x/dt = -k v v_x

dv_y/dt = -g - k v v_y

여기서 v = √(v_x² + v_y²). 레온하르트 오일러는 18세기에 이 시스템을 해석적으로 해결했습니다. 오늘날에는 룬게-쿠타 4–5차 수준의 수치 적분을 통해 접근합니다.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

class BallisticsSimulator:
    def __init__(self):
        self.g = 9.81
        self.rho = 1.225

    def calculate_k(self, diameter, mass, Cd=0.3):
        A = np.pi * (diameter / 2) ** 2
        k = (Cd * self.rho * A) / (2 * mass)
        return k

    def equations_of_motion(self, t, state, k):
        x, vx, y, vy = state
        v = np.sqrt(vx**2 + vy**2)
        ax = -k * v * vx
        ay = -self.g - k * v * vy
        return [vx, ax, vy, ay]

    def simulate(self, v0, angle_deg, x0=0, y0=0,
              diameter=0.1, mass=10, Cd=0.3,
              t_max=100, dt=0.01):
        k = self.calculate_k(diameter, mass, Cd)
        angle_rad = np.radians(angle_deg)
        vx0 = v0 * np.cos(angle_rad)
        vy0 = v0 * np.sin(angle_rad)
        initial_state = [x0, vx0, y0, vy0]
        t_span = (0, t_max)
        t_eval = np.arange(0, t_max, dt)
        solution = solve_ivp(
            self.equations_of_motion,
            t_span,
            initial_state,
            args=(k,),
            t_eval=t_eval,
            method='RK45'
        )
        t = solution.t
        x = solution.y[0]
        y = solution.y[2]
        fall_index = np.where(y < 0)[0]
        if len(fall_index) > 0:
            fall_index = fall_index[0]
            t = t[:fall_index + 1]
            x = x[:fall_index + 1]
            y = y[:fall_index + 1]
        return t, x, y, k

    def plot_trajectory(self, t, x, y, k, v0, angle):
        plt.figure(figsize=(12, 8))
        plt.plot(x, y, 'b-', linewidth=2, label='궤도')
        plt.scatter(x[0], y[0], color='red', s=100, zorder=5, label='발사')
        plt.scatter(x[-1], y[-1], color='green', s=100, zorder=5, label='충돌')
        plt.xlabel('거리, m', fontsize=12)
        plt.ylabel('고도, m', fontsize=12)
        plt.title(f'사격 궤도\nV₀ = {v0} m/s, 각도 = {angle}°, k = {k:.6f}',
                  fontsize=14)
        plt.grid(True, alpha=0.3)
        plt.legend()
        plt.axis('equal')
        plt.show()

    def print_results(self, t, x, y):
        print("=== 시뮬레이션 결과 ===")
        print(f"사거리: {x[-1]:.2f} m")
        print(f"최고 고도: {max(y):.2f} m")
        print(f"비행 시간: {t[-1]:.2f} s")

# 예시: v0=250 m/s, 각도=60°, 지름=0.2 m, 질량=20 kg

v₀ = 250 m/s, α = 60°, 지름 20 cm, 질량 20 kg, Cd = 0.3인 경우, 계수 k는 k = (Cd ρ A)/(2m)로 계산되며, 여기서 A는 단면적입니다. RK45 방법을 사용한 solve_ivp 함수는 충돌 지점(y = 0)까지 정확도를 보장합니다.

시뮬레이션 파라미터 및 결과

핵심 입력 파라미터:

• 초기 속도 v₀: 20–250 m/s
• 발사 각도 α: 30–60°
• 질량 m: 0.5–20 kg
• 지름: 10–20 cm
• Cd: 0.3 (표준 값)
• k: 자동 계산

제곱형 저항 사례의 결과:

• 비행 사거리: 약 10–20 km (파라미터에 따라 다름)
• 최고 고도: 5–15 km
• 비행 시간: 50–100초
• 최종 속도: 종단 속도 근처

각 사례를 비교하면 공기 저항을 고려할 경우 사거리가 30–50% 감소함을 알 수 있습니다.

핵심 요약

• 저항 없음 → 포물선 궤도; 최대 사거리는 α = 45°일 때 발생.
• 스토크스 법칙은 낮은 속도(<30 m/s)에 적용되며, 지수 감쇠 특성 있음.
• 뉴턴 법칙은 포병 분야에서 표준이며, 수치적 상미분 방정식 해법 필요.
• scipy.integrate를 활용한 파이썬 코드는 실제 Cd와 ρ 값을 쉽게 반영할 수 있어 유연함.

— Editorial Team