대한민국 공군·국방과학연구소 한국형 극초음속 활공비행체(HGV) 마하 7 활공 도약(Glide-Skip) 궤적 에너지 관리 및 공력열 제어 1,000Hz 0-Alloc C# 재돌입 비행제어 엔진

제1장: 3차원 구면 지구 대기권 재돌입 비행 동역학 방정식

지구 자전을 고려한 구면 좌표계에서 HGV 질점의 위치를 지심 반경 $r$, 지리학적 경도 $\lambda$, 위도 $\phi$, 대기 상대속도를 $V$, 비행경사각을 $\gamma$, 진행방위각을 $\psi$라 할 때의 3차원 재돌입 운동 방정식은 다음과 같습니다:

제2장: Sutton-Graves 공력열 유속 및 3중 비행 한계 회랑(Flight Corridor)

HGV가 마하 7로 비행할 때 기체 열보호시스템(TPS: Thermal Protection System)과 구조적 한계를 보호하기 위해 반드시 준수해야 하는 3대 물리적 구속 조건은 다음과 같습니다:

제3장: 1,000Hz 0-Alloc C# 극초음속 유도조종 엔진 아키텍처

1ms 실시간 루프에서 고도별 대기밀도 지수 모델, 공력계수, 열유속 한계 및 뱅크각 반전을 가비지 컬렉션(GC) 부하 없이 초고속 계산하는 순수 C# 아키텍처입니다:

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Aerospace.Hypersonic
{
    /// <summary>
    /// HGV 비행 상태 텔레메트리 (0-Alloc stack-only)
    /// </summary>
    public readonly struct HgvFlightState
    {
        public readonly float AltitudeM;       // 해수면 기준 고도 h (m)
        public readonly float VelocityMps;     // 대기 상대속도 V (m/s)
        public readonly float FlightPathRad;   // 비행경사각 gamma (rad)
        public readonly float HeadingRad;      // 진행방위각 psi (rad)
        public readonly float CrossTrackErrM;  // 기준 항법선 대비 횡방향 오차 (m)

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public HgvFlightState(float alt, float vel, float gamma, float psi, float crossErr)
        {
            AltitudeM = alt;
            VelocityMps = vel;
            FlightPathRad = gamma;
            HeadingRad = psi;
            CrossTrackErrM = crossErr;
        }
    }

    /// <summary>
    /// HGV 실시간 유도 명령 및 안전 감시 출력 구조체
    /// </summary>
    public readonly struct HgvGuidanceOutput
    {
        public readonly float CommandedAlphaRad; // 받음각(AoA) 명령 (rad)
        public readonly float CommandedBankRad;  // 뱅크각(Bank) 명령 (rad)
        public readonly float HeatFluxKwM2;      // 정체점 열유속 (kW/m^2)
        public readonly float DynamicPressureKpa;// 동압 (kPa)
        public readonly float LoadFactorG;       // 구조 하중 (G)
        public readonly bool IsCorridorSafe;     // 3대 비행 한계 회랑 만족 여부

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public HgvGuidanceOutput(float alpha, float bank, float qHeat, float qDyn, float nz, bool safe)
        {
            CommandedAlphaRad = alpha;
            CommandedBankRad = bank;
            HeatFluxKwM2 = qHeat;
            DynamicPressureKpa = qDyn;
            LoadFactorG = nz;
            IsCorridorSafe = safe;
        }
    }

    /// <summary>
    /// 공군·ADD 극초음속 활공비행체(HGV) 1,000Hz 0-Alloc 재돌입 유도조종 코어
    /// </summary>
    public ref struct HgvReentryGuidanceEngine
    {
        // 물리 및 기체 제원 상수
        public const float EarthRadiusM = 6371000.0f;
        public const float SeaLevelGravity = 9.80665f;
        public const float Rho0 = 1.225f;           // 해수면 공기밀도 (kg/m^3)
        public const float ScaleHeightM = 7200.0f;   // 대기 밀도 스케일 높이 H_s (m)
        public const float MassKg = 1200.0f;        // HGV 비행체 질량 (1.2톤)
        public const float RefAreaM2 = 0.85f;       // 기준 단면적 S_ref (m^2)
        public const float NoseRadiusM = 0.05f;     // 노즈 노즐 반경 R_N (5cm)
        public const float Ksg = 1.7415e-4f;        // Sutton-Graves 공력열 계수
        public const float Vco = 7900.0f;           // 제1우주원궤도속도 (m/s)

        // 안전 제약 한계치 (Corridor Limits)
        public const float MaxHeatFluxKwM2 = 2500.0f;
        public const float MaxDynamicPressureKpa = 150.0f;
        public const float MaxLoadFactorG = 6.0f;
        public const float CrossTrackCorridorWidthM = 25000.0f; // 횡방향 회랑 25km

        /// <summary>
        /// 1,000Hz 결정론적 극초음속 준평형 활공 유도(QEGC) 및 뱅크각 반전 명령 산출
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static HgvGuidanceOutput EvaluateGuidance(in HgvFlightState state, float currentBankSign)
        {
            float r = EarthRadiusM + state.AltitudeM;
            float g = SeaLevelGravity * (EarthRadiusM * EarthRadiusM) / (r * r);

            // 1. 고도별 대기 밀도 지수 모델 계산
            float rho = Rho0 * MathF.Exp(-state.AltitudeM / ScaleHeightM);

            // 2. 동압(Dynamic Pressure) 및 정체점 열유속(Sutton-Graves) 계산
            float v = state.VelocityMps;
            float qDynPa = 0.5f * rho * v * v;
            float qDynKpa = qDynPa * 0.001f;

            float normRho = MathF.Max(0.000001f, rho / Rho0);
            float normV = v / Vco;
            float heatFluxKwM2 = (Ksg / MathF.Sqrt(NoseRadiusM)) * MathF.Sqrt(normRho) * MathF.Pow(normV, 3.15f) * 1000.0f;

            // 3. 받음각(Alpha) 프로파일 결정 (속도에 따른 양항비 최적화)
            // 고속(마하 7, 2,400m/s)에서는 고받음각(18도)으로 에너지 감속, 저속에서는 최대 양항비(10도) 활공
            float alphaRad = (v > 1800.0f) ? 0.314f : 0.174f; // 18 deg -> 10 deg

            // 극초음속 수정 뉴턴 공력 계수
            float cl = 0.10f + 1.20f * alphaRad;
            float cd = 0.05f + 1.60f * alphaRad * alphaRad;
            float liftN = qDynPa * RefAreaM2 * cl;
            float dragN = qDynPa * RefAreaM2 * cd;

            float loadFactorG = MathF.Sqrt(liftN * liftN + dragN * dragN) / (MassKg * SeaLevelGravity);

            // 4. 준평형 활공 유도(QEGC) 뱅크각 크기 계산
            // L * cos(sigma) = m * (g - v^2 / r)
            float requiredVerticalForce = MassKg * (g - (v * v / r));
            float cosBank = Math.Clamp(requiredVerticalForce / MathF.Max(1.0f, liftN), 0.0f, 1.0f);
            float bankMagRad = MathF.Acos(cosBank);

            // 열유속 초과 임박 시 기수를 들어 대기권 상층으로 궤적 완화 (동적 풀업)
            if (heatFluxKwM2 > MaxHeatFluxKwM2 * 0.90f)
            {
                bankMagRad *= 0.60f; // 뱅크각을 줄여 수직 양력을 확보함으로써 고도 하강 억제
            }

            // 5. 횡방향 에너지 관리를 위한 뱅크각 반전(Bank-Angle Reversal) 로직
            float commandedBankSign = currentBankSign;
            if (state.CrossTrackErrM > CrossTrackCorridorWidthM)
            {
                commandedBankSign = -1.0f; // 우측 이탈 시 좌측 선회로 반전
            }
            else if (state.CrossTrackErrM < -CrossTrackCorridorWidthM)
            {
                commandedBankSign = 1.0f;  // 좌측 이탈 시 우측 선회로 반전
            }

            float finalBankCmdRad = bankMagRad * commandedBankSign;

            // 6. 3대 안전 회랑 종합 판정
            bool isCorridorSafe = (heatFluxKwM2 <= MaxHeatFluxKwM2) &&
                                 (qDynKpa <= MaxDynamicPressureKpa) &&
                                 (loadFactorG <= MaxLoadFactorG);

            return new HgvGuidanceOutput(
                alphaRad,
                finalBankCmdRad,
                heatFluxKwM2,
                qDynKpa,
                loadFactorG,
                isCorridorSafe
            );
        }
    }
}

제4장: 2022 개정 교육과정 연계 교과 융합 지도안 (교사용 가이드)