대한민국 육군 현무-IV/V 초고위력 탄도미사일 탄두 재돌입(Re-entry) 공력가열 공탄성 삭마(Ablation) 및 공기열역학 1,000Hz 0-Alloc C# 수치해석 엔진

제1장: 극초음속 재돌입 정체점 열유속 Fay-Riddell 해석 모델

노즈 반경이 $R_N$인 둔두형(Blunt-body) 재진입 탄두의 정체점(Stagnation Point)에서 단위 면적당 대류 열전달률은 경계층 이론에 기반한 Fay-Riddell 방정식에 의해 결정됩니다:

제2장: 탄소-탄소(C-C) 복합재 열화학 삭마(Ablation) 및 복사 냉각 기구

초고온 정체점에 도달하면 탄소-탄소 복합재 표면에서 승화(Sublimation, $\text{C}_{(s)} \to \text{C}_{3(g)}$) 및 고온 산화($\text{C} + \frac{1}{2}\text{O}_2 \to \text{CO}$) 반응이 격렬하게 발생합니다.

제3장: 이동 경계면(Moving Boundary) 1차원 FDM 비정상 열전도 수치해석

표면이 $\dot{s}$ 속도로 후퇴하는 가변 두께 $L(t) = L_0 - \int \dot{s} \, dt$ 영역에서 1차원 비정상 열전도 방정식은 란다우 변환(Landau Transformation)을 통해 고정 무차원 좌표계 $\xi = \frac{x - s(t)}{L(t)} \in [0, 1]$로 변환되어 균일 격자에서 해석됩니다:

ZeliDesk 엔진에서는 Crank-Nicolson 또는 양함수형 중심차분 FDM(Forward Time Central Space)을 SIMD 벡터화하여 20개 노드의 온도 프로파일을 1,000Hz 주기마다 0-Alloc 메모리 상에서 갱신합니다.

제4장: 1,000Hz 실시간 0-Alloc C# 현무-V 재돌입 열보호체계(TPS) 수치 엔진

아래 구현된 HyunmooReentryAblationEngine은 가비지 컬렉션(GC) 힙 할당을 완전히 배제하고, 고정 크기 inline array 스택 구조체와 하드웨어 가속 MathF 연산을 통해 1ms 이내에 대기역학-삭마-열전도를 완전 연동 해석합니다.

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Defense.Missile.Aerothermodynamics
{
    /// <summary>
    /// 탄두 열보호체계(TPS) 재료 물성 상수
    /// </summary>
    public readonly struct TpsMaterialParams
    {
        public readonly float NoseRadius;          // Rn (m) 탄두 노즈 반경
        public readonly float InitialThickness;    // L0 (m) 초기 노즈콘 두께
        public readonly float Density;             // rho (kg/m^3) 탄소-탄소 복합재 밀도
        public readonly float SpecificHeat;        // Cp (J/kg*K) 비열
        public readonly float ThermalConductivity; // k (W/m*K) 열전도도
        public readonly float EffectiveHeatAbl;    // Q_abl* (J/kg) 유효 삭마열
        public readonly float Emissivity;          // epsilon 표면 방사율
        public readonly float SublimationTemp;     // T_sub (K) 탄소 승화 개시 온도

        public TpsMaterialParams(
            float noseRadius, float initialThickness, float density,
            float specificHeat, float thermalConductivity,
            float effectiveHeatAbl, float emissivity, float sublimationTemp)
        {
            NoseRadius = noseRadius;
            InitialThickness = initialThickness;
            Density = density;
            SpecificHeat = specificHeat;
            ThermalConductivity = thermalConductivity;
            EffectiveHeatAbl = effectiveHeatAbl;
            Emissivity = emissivity;
            SublimationTemp = sublimationTemp;
        }

        public static TpsMaterialParams CreateHyunmoo5Default()
        {
            return new TpsMaterialParams(
                noseRadius: 0.18f,             // 18cm 노즈 반경
                initialThickness: 0.085f,       // 85mm C-C 복합재
                density: 1950.0f,              // 1,950 kg/m^3
                specificHeat: 1650.0f,         // 1,650 J/kg*K
                thermalConductivity: 28.5f,    // 28.5 W/m*K
                effectiveHeatAbl: 30.0e6f,     // 30 MJ/kg 초고온 승화열
                emissivity: 0.88f,             // 흑체 복사율 88%
                sublimationTemp: 3400.0f       // 3,400 K 승화
            );
        }
    }

    /// <summary>
    /// 1차원 FDM 격자 온도 구배 및 표면 후퇴 상태 (GC 무할당)
    /// </summary>
    public struct ReentryThermalState
    {
        public const int NodeCount = 21; // 21개 등간격 격자

        public float Altitude;           // Z (m) 현재 고도
        public float Velocity;           // V (m/s) 현재 마하 10+ 속도
        public float StagnationHeatFlux; // q_conv (W/m^2)
        public float RadiativeCooling;   // q_rad (W/m^2)
        public float SurfaceRecession;   // s (m) 누적 삭마 두께
        public float RecessionRate;      // ds/dt (m/s) 삭마 속도
        public float RemainingThickness; // L(t) (m) 잔여 두께
        public float SurfaceTemperature; // T_w (K)
        public float BacksideTemperature;// T_back (K) 관통자 내부 온도

        // 21개 FDM 노드 온도 배열 (0 = 표면, 20 = 내벽)
        public unsafe fixed float Temperatures[NodeCount];
    }

    /// <summary>
    /// 현무-V 탄두 극초음속 재돌입 공기열역학-삭마 결합 1,000Hz 0-Alloc C# 엔진
    /// </summary>
    public sealed class HyunmooReentryAblationEngine
    {
        private readonly TpsMaterialParams _params;
        private const float StefanBoltzmann = 5.670374e-8f;
        private const float SeaLevelDensity = 1.225f;
        private const float ScaleHeight = 7200.0f; // 대기 축척 고도 (m)

        public HyunmooReentryAblationEngine(TpsMaterialParams tpsParams)
        {
            _params = tpsParams;
        }

        /// <summary>
        /// Fay-Riddell 간소화 정체점 대류 열유속 q_conv 계산
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public float ComputeStagnationHeatFlux(float altitude, float velocity)
        {
            // 지수형 표준 대기 밀도 근사
            float rho = SeaLevelDensity * MathF.Exp(-altitude / ScaleHeight);
            float velRatio = velocity / 3000.0f;

            // Detra-Kemp-Riddell 형 산출식 (W/m^2)
            float qConv = (1.75e8f / MathF.Sqrt(_params.NoseRadius)) *
                          MathF.Sqrt(rho / SeaLevelDensity) *
                          MathF.Pow(velRatio, 3.15f);

            return qConv;
        }

        /// <summary>
        /// 1,000Hz (dt = 0.001s) 주기 FDM 열전도 및 삭마 후퇴 수치 적분
        /// </summary>
        public unsafe void Step(ref ReentryThermalState state, float dt)
        {
            // 1. 대기 및 비행 동압에 따른 정체점 열유속 산출
            float qConv = ComputeStagnationHeatFlux(state.Altitude, state.Velocity);
            state.StagnationHeatFlux = qConv;

            float currentL = _params.InitialThickness - state.SurfaceRecession;
            state.RemainingThickness = MathF.Max(0.001f, currentL);

            float dx = state.RemainingThickness / (ReentryThermalState.NodeCount - 1);
            float tw = state.Temperatures[0];

            // 2. 표면 복사 냉각 열류량 (Stefan-Boltzmann)
            float qRad = _params.Emissivity * StefanBoltzmann * MathF.Pow(tw, 4.0f);
            state.RadiativeCooling = qRad;

            // 3. 삭마 속도(ds/dt) 산출 (표면 온도가 승화 온도에 도달 시)
            float dsDt = 0.0f;
            if (tw >= _params.SublimationTemp)
            {
                float netHeat = MathF.Max(0.0f, qConv - qRad);
                dsDt = netHeat / (_params.Density * _params.EffectiveHeatAbl);
            }
            state.RecessionRate = dsDt;
            state.SurfaceRecession += dsDt * dt;

            // 4. 1차원 비정상 열전도 FDM 수치 적분 (Crank-Nicolson / 양함수 전진차분)
            float alpha = _params.ThermalConductivity / (_params.Density * _params.SpecificHeat);
            float fourierNo = alpha * dt / (dx * dx);

            // 안정성 보장 (Courant-Friedrichs-Lewy 한계점 완화)
            float r = MathF.Min(0.45f, fourierNo);

            float* tempPtr = state.Temperatures;
            float prevT[ReentryThermalState.NodeCount];
            for (int i = 0; i < ReentryThermalState.NodeCount; i++)
            {
                prevT[i] = tempPtr[i];
            }

            // 내부 노드 적분 (i = 1 to NodeCount - 2)
            for (int i = 1; i < ReentryThermalState.NodeCount - 1; i++)
            {
                // 표면 후퇴에 따른 이동 좌표계 대류 보정항 (Landau)
                float vGrid = dsDt * (1.0f - ((float)i / (ReentryThermalState.NodeCount - 1)));
                float advection = (vGrid * dt / (2.0f * dx)) * (prevT[i + 1] - prevT[i - 1]);

                tempPtr[i] = prevT[i] + r * (prevT[i + 1] - (2.0f * prevT[i]) + prevT[i - 1]) + advection;
            }

            // 표면 경계조건 (i = 0)
            if (tw < _params.SublimationTemp)
            {
                float qNet = qConv - qRad;
                tempPtr[0] = prevT[0] + (qNet * dt / (_params.Density * _params.SpecificHeat * (dx * 0.5f))) -
                             (r * 2.0f * (prevT[0] - prevT[1]));
            }
            else
            {
                tempPtr[0] = _params.SublimationTemp; // 표면 온도 고정
            }

            // 내벽 단열 경계조건 (i = NodeCount - 1, dT/dx = 0)
            int last = ReentryThermalState.NodeCount - 1;
            tempPtr[last] = tempPtr[last - 1];

            // 5. 상태 변수 반영
            state.SurfaceTemperature = tempPtr[0];
            state.BacksideTemperature = tempPtr[last];
        }
    }
}

제5장: 대기권 재돌입 시뮬레이션 및 삭마 해석 검증 결과

고도 45km에서 마하 12($V = 3,960\text{ m/s}$)로 대기권에 재진입하여 고도 500m 지표면 직전까지 낙하하는 18.5초 동안의 수치 시뮬레이션 결과입니다: