Python에서 Cx(M)와 ρ(h)를 사용한 탄도 모델링

가변 공기 밀도를 고려한 외부 탄도학 수치 모델링

외부 탄도학 문제 분석을 이어가며, 항력 계수 Cx가 속도에 의존하고 공기 밀도가 고도에 따라 변하는 경우를 살펴봅니다. 네 번째와 다섯 번째 경우는 해석적 해가 불가능하므로 수치적 방법이 필요합니다. 발사체의 비행은 초기 속도 600 m/s, 각도 55°, 질량 50 kg의 조건으로 시뮬레이션됩니다.

1943년 항력 법칙은 Cx를 마하수 M = v/a의 함수로 정의합니다. 밀도 ρ(h)는 지수 함수 ρ(h) = ρ₀ × exp(-h/H)로 근사하며, 여기서 H = 8000 m는 척도 높이입니다.

마하수에 따른 Cx 의존성

계수 Cx는 구간별 함수입니다:

• M ≤ 0.73: Cx = 0.157
• 0.73 < M < 0.82: Cx = 0.033M + 0.133
• 0.82 ≤ M < 0.91: Cx = 0.161 + 3.9(M - 0.823)²
• 0.91 ≤ M < 1.00: Cx = 1.5M - 1.176
• M ≥ 1.00: Cx = 0.384 - 1.6(M - 1.176)²

항력 F = 0.5 ρ(h) v² S Cx (v/|v|)이며, S는 단면적(0.01 m²)입니다. 운동 방정식:

ẋ = vx

ẏ = vy

v̇x = -Fx/m

v̇y = -g - Fy/m

g = 9.81 m/s², a = 343 m/s.

수치 적분 구현

상미분 방정식 시스템을 풀기 위해 4차 룽게-쿠타 방법(RK4)을 적용합니다. 적분 간격 dt = 0.01 s, 최종 시간 t_end = 98.4 s.

Python 전체 시뮬레이션 코드:

import numpy as np
import matplotlib.pyplot as plt

# 상수
g = 9.81  # m/s²
m = 50.0  # 질량, kg
S = 0.01  # 단면적, m²
rho0 = 1.225  # 해수면 공기 밀도, kg/m³
H = 8000  # 척도 높이, m
a = 343  # 음속, m/s
v0 = 600  # 초기 속도, m/s
theta = np.radians(55)  # 초기 각도, rad
dt = 0.01  # 적분 간격, s
t_end = 98.4  # 최종 시간, s

# 초기 조건
x, y = 0, 0
vx, vy = v0 * np.cos(theta), v0 * np.sin(theta)

# 결과 저장 배열
t = np.arange(0, t_end, dt)
X, Y = np.zeros_like(t), np.zeros_like(t)
VX, VY = np.zeros_like(t), np.zeros_like(t)


# cx(M) 계산 함수
def cx(M):
    if M <= 0.73:
        return 0.157
    elif M < 0.82:
        return 0.033 * M + 0.133
    elif M < 0.91:
        return 0.161 + 3.9 * (M - 0.823) ** 2
    elif M < 1.00:
        return 1.5 * M - 1.176
    else:
        return 0.384 - 1.6 * (M - 1.176) ** 2


# 공기 밀도 함수
def rho(h):
    return rho0 * np.exp(-h / H)


# 상미분 방정식 시스템
def f(state, t):
    x, y, vx, vy = state
    v = np.sqrt(vx ** 2 + vy ** 2)
    M = v / a
    rho_h = rho(y)
    cx_val = cx(M)

    Fx = 0.5 * rho_h * v ** 2 * S * cx_val * vx / v
    Fy = 0.5 * rho_h * v ** 2 * S * cx_val * vy / v

    dxdt = vx
    dydt = vy
    dvxdt = -Fx / m
    dvydt = -g - Fy / m

    return np.array([dxdt, dydt, dvxdt, dvydt])


# 4차 룽게-쿠타 방법
def rk4_step(state, t, dt):
    k1 = dt * f(state, t)
    k2 = dt * f(state + 0.5 * k1, t + 0.5 * dt)
    k3 = dt * f(state + 0.5 * k2, t + 0.5 * dt)
    k4 = dt * f(state + k3, t + dt)
    return state + (k1 + 2 * k2 + 2 * k3 + k4) / 6


# 적분
state = np.array([x, y, vx, vy])
for i, t_val in enumerate(t):
    X[i], Y[i], VX[i], VY[i] = state
    state = rk4_step(state, t_val, dt)

# 시각화
plt.figure(figsize=(10, 6))
plt.plot(X, Y, label='비행 궤적')
plt.xlabel('X, m')
plt.ylabel('Y, m')
plt.title('발사체 비행 궤적')
plt.grid(True)
plt.legend()
plt.show()

# 결과 출력
print(f"비행 거리: {X[-1]:.2f} m")
print(f"최대 고도: {np.max(Y):.2f} m")

시뮬레이션 결과

비행 거리: 29 km. 최대 고도: 11.5 km. 궤적은 초음속 속도 전환과 고도에 따른 공기 희박화를 고려합니다.

주요 가정:

• 지수 밀도 모델(단순화됨; 실제 대기에는 더 정확한 프로필 필요).
• 다른 효과(회전, 바람)를 고려하지 않은 2차 항력.

| 매개변수 | 값 |

|-----------|-------|

| 비행 거리 | 29,000 m |

| 최대 고도 | 11,500 m |

| 비행 시간 | 98.4 s |

| 초기 속도 | 600 m/s |

핵심 요약

• 수치적 방법은 필수적입니다: 가변 Cx(M)와 ρ(h)를 포함한 현실적 모델링에 필요.
• RK4는 정확도 보장: 탄도학 작업에서 dt=0.01 s로 충분한 정밀도 제공.
• 1943년 법칙 정확도: 천음속 영역(M≈0.8–1.0)을 정확히 모델링.
• 지수 ρ(h) 합리적 근사: 15 km까지의 고도에서 타당한 근사치 제공.

— Editorial Team