6DOF 로켓 시뮬레이션

6DOF란 무엇인가?

6DOF는 6 Degrees of Freedom의 약자로, 한국어로는 "6자유도"라고 합니다. 3차원 공간에서 물체가 가질 수 있는 모든 운동의 종류를 의미합니다.

병진 운동 (Translation) - 3 DOF

물체의 위치가 변하는 운동입니다. X, Y, Z 세 방향으로 이동할 수 있습니다.

회전 운동 (Rotation) - 3 DOF

물체의 자세가 변하는 운동입니다. 세 축을 중심으로 회전할 수 있습니다.

왜 6DOF 시뮬레이션이 필요한가?

간단한 시뮬레이션에서는 로켓이 수직으로만 올라간다고 가정합니다 (1DOF). 하지만 실제 로켓은:

• 바람에 의해 옆으로 밀립니다
• 중심이동(CG)과 압력중심(CP)의 차이로 회전합니다
• 추력이 완벽히 중심을 통과하지 않아 비틀림이 생깁니다
• 공기저항이 방향에 따라 다르게 작용합니다

이러한 복잡한 현상을 정확히 예측하려면 6DOF 시뮬레이션이 필수입니다.

좌표계 시스템

로켓 시뮬레이션에서는 여러 좌표계를 사용합니다. 각 좌표계는 서로 다른 목적에 최적화되어 있습니다.

NED 좌표계 (관성 좌표계)

% 상태벡터에서 위치 추출
pos = X(1:3);  % [X_north; Y_east; Z_down]

% 고도 계산 (NED에서 Z는 아래가 +)
alt = -pos(3);  % 고도 = -Z_down

Body 좌표계 (기체 좌표계)

Body frame은 공력 계산에 유용합니다. 로켓 입장에서 바람이 어느 방향에서 오는지 쉽게 알 수 있기 때문입니다.

좌표계 변환

두 좌표계 사이의 변환은 회전 행렬(DCM, Direction Cosine Matrix)을 사용합니다.

여기서:

• \(\vec{v}_{NED}\): NED 좌표계에서의 벡터
• \(\vec{v}_{Body}\): Body 좌표계에서의 벡터
• \(\mathbf{R}_{B \to I}\): Body에서 관성(Inertial) 좌표계로의 회전 행렬

% Body frame 속도를 NED frame으로 변환
R_B2I = GetDCM_QUAT(quat);  % 쿼터니언에서 회전행렬 계산
V_NED = R_B2I * V_Body;    % 좌표 변환

% 중력을 Body frame으로 변환
R_I2B = R_B2I';            % 역변환 = 전치행렬
g_I = [0; 0; 9.81];      % NED에서 중력 (아래로 +)
F_grav = R_I2B * (m * g_I); % Body frame 중력

쿼터니언과 자세 표현

3D 회전을 표현하는 방법은 여러 가지가 있습니다. 그 중 쿼터니언(Quaternion)이 시뮬레이션에서 가장 많이 사용됩니다.

자세 표현 방법 비교

짐벌락 (Gimbal Lock) 문제

쿼터니언의 정의

쿼터니언은 4개의 요소로 구성된 수입니다:

여기서:

• \(q_w\): 스칼라(실수) 부분
• \(q_x, q_y, q_z\): 벡터(허수) 부분
• \(\mathbf{i}, \mathbf{j}, \mathbf{k}\): 허수 단위 (i² = j² = k² = ijk = -1)

단위 쿼터니언

회전을 표현하려면 단위 쿼터니언을 사용해야 합니다:

% 시뮬레이션 중 쿼터니언 정규화 (크기를 1로 유지)
X(7:10) = X(7:10) / norm(X(7:10));

오일러 각에서 쿼터니언으로 변환

ZYX 순서 (Yaw → Pitch → Roll)의 오일러 각을 쿼터니언으로 변환하는 공식:

function quat = GetQUAT(psi, theta, phi)
    % 삼각함수 미리 계산 (효율성)
    cPhi = cos(phi/2); sPhi = sin(phi/2);
    cTheta = cos(theta/2); sTheta = sin(theta/2);
    cPsi = cos(psi/2); sPsi = sin(psi/2);
    
    % 쿼터니언 계산
    qw = cPhi*cTheta*cPsi + sPhi*sTheta*sPsi;
    qx = sPhi*cTheta*cPsi - cPhi*sTheta*sPsi;
    qy = cPhi*sTheta*cPsi + sPhi*cTheta*sPsi;
    qz = cPhi*cTheta*sPsi - sPhi*sTheta*cPsi;
    
    quat = [qw, qx, qy, qz];
end

쿼터니언에서 회전 행렬로 변환

쿼터니언 미분 (시간에 따른 변화)

각속도 \(\vec{\omega} = [\omega_x, \omega_y, \omega_z]^T\)가 주어졌을 때, 쿼터니언의 시간 미분은:

여기서 \(\Omega\)는 스큐-대칭 행렬(Skew-symmetric matrix)입니다:

function dot_q = Derivative_Quat(q, omega_b)
    q = q(:);          % 열벡터로 변환
    omega_b = omega_b(:);
    
    % Omega 행렬 구성
    Omega = [    0,      -omega_b(1), -omega_b(2), -omega_b(3);
             omega_b(1),       0,      omega_b(3), -omega_b(2);
             omega_b(2), -omega_b(3),        0,     omega_b(1);
             omega_b(3),  omega_b(2), -omega_b(1),        0   ];
    
    dot_q = 0.5 * Omega * q;
end

뉴턴 역학과 운동방정식

로켓의 운동은 뉴턴의 운동 법칙을 따릅니다. 병진 운동과 회전 운동 각각에 대한 방정식이 필요합니다.

뉴턴의 제2법칙 (병진 운동)

Body frame에서 속도 미분을 계산할 때는 회전 효과를 고려해야 합니다:

여기서:

• \(\vec{v}_B\): Body frame에서의 속도
• \(\vec{F}_{total}\): 총 외력 (추력 + 공력 + 중력)
• \(\vec{\omega}\): 각속도 벡터
• \(\vec{\omega} \times \vec{v}_B\): 코리올리 효과 (회전하는 좌표계의 효과)

오일러 방정식 (회전 운동)

회전 운동에 대한 뉴턴 제2법칙입니다:

여기서:

• \(\vec{M}\): 총 모멘트 (토크)
• \(\mathbf{J}\): 관성 모멘트 텐서 (3×3 행렬)
• \(\dot{\vec{\omega}}\): 각가속도

각가속도를 구하면:

관성 모멘트 텐서

로켓을 원통으로 근사하면, 관성 모멘트는 대각 행렬이 됩니다:

원통의 관성 모멘트:

% 관성모멘트 (원통 근사)
rocket.Ixx = 0.5 * m0 * (rocket.d/2)^2;            % 롤 축
rocket.Iyy = m0 * (rocket.L^2/12 + (rocket.d/2)^2/4); % 피치 축
rocket.Izz = rocket.Iyy;                              % 요 축 (대칭)
rocket.J = diag([rocket.Ixx, rocket.Iyy, rocket.Izz]);

로켓에 작용하는 힘들

공기역학 (Aerodynamics)

로켓이 대기 중을 날아갈 때, 공기와의 상호작용으로 여러 힘과 모멘트가 발생합니다.

받음각 (Angle of Attack, α)

옆미끄럼각(Sideslip angle, β)도 비슷하게 계산합니다:

% 받음각 계산 (Body frame 기준)
if V > 0.1
    alpha = atan2(V_B(3), V_B(1));  % 피치 평면
    beta = atan2(V_B(2), V_B(1));   % 요 평면
else
    alpha = 0;
    beta = 0;
end

동압 (Dynamic Pressure)

공기역학적 힘의 크기를 결정하는 중요한 량입니다:

여기서:

• \(\rho\): 공기 밀도 [kg/m³]
• \(V\): 속력 [m/s]

항력 (Drag)

여기서:

• \(A_{ref}\): 기준 면적 (로켓 단면적) = \(\pi r^2\)
• \(C_D\): 항력 계수

항력 계수는 받음각에 따라 증가합니다:

% 동압 계산
q_dyn = 0.5 * rho * V^2;

% 항력 계수 (받음각 의존성)
Cd = rocket.Cd0 * (1 + 0.3 * alpha^2);

% 항력 크기
D = q_dyn * rocket.A_ref * Cd;

% 항력 벡터 (속도 반대 방향)
F_aero = -D * (V_B / V);

법선력과 양력

받음각이 생기면 로켓 측면에 힘이 작용합니다:

여기서 \(C_{N\alpha}\)는 법선력 기울기(Normal force slope)로, 받음각 1 rad당 법선력 계수 변화량입니다.

마하수

속도를 음속으로 나눈 무차원 수입니다:

여기서 \(a\)는 음속으로, 고도에 따라 변합니다.

추진 시스템

로켓 모터는 추진제를 연소시켜 추력을 발생시킵니다. 이 시뮬레이션에서는 KTM I386 모터를 사용합니다.

추력 곡선

추력은 시간에 따라 변합니다. 모터 데이터로부터 추력 곡선을 정의합니다:

% 시간 [초]
motor.t = [0.0, 0.08, 0.16, ... , 2.30];

% 추력 [N]
motor.thrust = [0.0, 2.09, 19.78, ... , 0.0];

% 연소 시간
motor.t_burn = 2.30;  % [초]

% 현재 시간의 추력 보간
T_current = interp1(motor.t, motor.thrust, tc, 'linear', 0);

총 임펄스

총 임펄스는 추력을 시간에 대해 적분한 값입니다:

이 값은 모터의 "등급"을 결정합니다 (I급 모터 = 320~640 Ns).

% 수치 적분 (사다리꼴 법칙)
total_impulse = trapz(motor.t, motor.thrust);
fprintf('총 임펄스: %.1f Ns\n', total_impulse);

질량 변화

연소 중에 추진제가 소모되어 질량이 감소합니다:

if tc < motor.t_burn
    % 연소 중: 추진제 선형 감소
    prop_remaining = rocket.m_prop * (1 - tc/motor.t_burn);
    m_current = rocket.m_dry + rocket.m_case + prop_remaining;
    
    % CG 이동 (추진제 소모에 따라 앞으로)
    x_cg = rocket.x_cg0 + 0.05 * (tc/motor.t_burn);
else
    % 연소 종료 후: 고정
    m_current = rocket.m_dry + rocket.m_case;
    x_cg = rocket.x_cg0 + 0.05;
end

대기 모델

공기의 밀도, 온도, 음속은 고도에 따라 변합니다. MATLAB의 atmosisa 함수는 국제 표준 대기(ISA)를 사용합니다.

국제 표준 대기 (ISA)

압력과 밀도

음속

이상기체에서 음속은:

여기서:

• \(\gamma = 1.4\) (공기의 비열비)
• \(R = 287\) J/(kg·K) (공기의 기체상수)

% 현재 고도의 대기 속성
alt = -pos(3);  % NED에서 고도
[~, a, ~, rho] = atmosisa(max(0, alt));

% a: 음속 [m/s]
% rho: 공기 밀도 [kg/m³]

% 마하수 계산
Mach = V / a;

안정성과 댐핑

로켓이 안정적으로 비행하려면 자세가 흐트러져도 다시 원래대로 돌아와야 합니다.

안정 마진 (Static Margin)

복원 모멘트

안정한 로켓이 받음각을 가지면, 원래 방향으로 되돌리려는 모멘트가 발생합니다:

음수 부호는 받음각을 줄이는 방향으로 작용함을 의미합니다.

% 안정 마진
SM = rocket.x_cp - x_cg;  % 양수 = 안정

% 법선력
CN = rocket.CN_alpha * alpha;
F_N = q_dyn * rocket.A_ref * CN;

% 복원 모멘트 (피치)
M_restore_y = -F_N * SM;

% 복원 모멘트 (요)
M_restore_z = -q_dyn * rocket.A_ref * rocket.CN_alpha * beta * SM;

댐핑 모멘트

회전하는 로켓은 공기저항에 의해 회전이 감소합니다. 이것이 댐핑 효과입니다:

여기서:

• \(C_{mq}\): 피치 댐핑 계수 (보통 음수, 예: -10)
• \(d\): 로켓 직경
• \(\omega\): 각속도

if V > 1
    damp_factor = rocket.Cmq * q_dyn * rocket.A_ref * rocket.d^2 / (2*V);
    M_damp = [0; ...
              damp_factor * omega(2); ...  % 피치
              damp_factor * omega(3)];     % 요
else
    M_damp = [0; 0; 0];
end

수치 적분 (RK4)

미분방정식을 컴퓨터로 풀려면 수치적 방법이 필요합니다. 4차 Runge-Kutta (RK4) 방법이 가장 널리 사용됩니다.

오일러 방법 vs RK4

RK4 알고리즘

미분방정식 \(\frac{dy}{dt} = f(t, y)\)에 대해:

% RK4 적분
[k1, aero] = rocket_dynamics(X, T_current, m_current, x_cg, rocket);
[k2, ~] = rocket_dynamics(X + 0.5*dt*k1, T_current, m_current, x_cg, rocket);
[k3, ~] = rocket_dynamics(X + 0.5*dt*k2, T_current, m_current, x_cg, rocket);
[k4, ~] = rocket_dynamics(X + dt*k3, T_current, m_current, x_cg, rocket);

% 상태 업데이트
X = X + (dt/6) * (k1 + 2*k2 + 2*k3 + k4);

% 쿼터니언 정규화 (크기 1 유지)
X(7:10) = X(7:10) / norm(X(7:10));

시간 스텝 선택

시뮬레이션 주파수 1000 Hz (dt = 0.001 s)를 사용합니다. 이는 로켓 동역학의 시간 스케일에 비해 충분히 작습니다.

상태 벡터 구성

시뮬레이션에서 추적해야 할 모든 변수를 하나의 벡터로 모읍니다. 이것이 상태 벡터(State Vector)입니다.

13차원 상태 벡터

상태 미분 벡터

동역학 함수는 상태의 시간 미분을 반환합니다:

% 상태 미분
pos_dot = R_B2I * V_B;                              % 위치 미분
vel_dot = F_total/m - cross(omega, V_B);            % 속도 미분
quat_dot = Derivative_Quat(quat, omega);            % 쿼터니언 미분
omega_dot = J \ (M_total - cross(omega, J*omega));  % 각속도 미분

% 상태 미분 벡터
Xdot = [pos_dot; vel_dot; quat_dot; omega_dot];

초기 조건

% 초기 위치 (발사대)
pos0 = [0; 0; 0];

% 초기 속도 (정지)
vel0 = [0; 0; 0];

% 발사 각도 (86도 = 수직에서 4도 기울임)
launch_angle = 86;
euler0 = [0; deg2rad(launch_angle); 0];

% 오일러 각을 쿼터니언으로 변환
quat0 = GetQUAT(euler0(3), euler0(2), euler0(1))';

% 초기 각속도 (정지)
omega0 = [0; 0; 0];

% 전체 초기 상태 벡터
X0 = [pos0; vel0; quat0; omega0];

핵심 요약