대한민국 육군 극초음속 전술탄도탄 탄두 첨두부 탄소-탄소(C/C) 복합재 삭마(Ablation) 열전도 및 1,000Hz 0-Alloc C# 공력 형상 변화 보상 비행제어 엔진

1. 극초음속 재진입 공력 가열 및 정체점 열유속 역학

극초음속 비행체 전방의 뭉툭한 첨두부 주위에는 강력한 이탈 충격파(Detached Bow Shock)가 형성됩니다. 충격파 직후의 정체점(Stagnation Point)에서 자유류 기동 에너지가 열에너지로 변환되며, 화학적 비평형(Chemical Non-equilibrium) 해리 반응(Dissociation)이 발생합니다.

정체점 표면에서의 대류 열유속 $q_s$는 고전적인 Fay-Riddell 공식에 의해 다음과 같이 기술됩니다:

여기서 $\mathrm{Pr}$은 프란틀 수, $ ho_w, \mu_w$는 표면 벽면에서의 기체 밀도 및 점성계수, $ ho_e, \mu_e$는 경계층 외곽에서의 밀도 및 점성계수, $(du_e/dx)_s$는 정체점 속도 구배(Velocity Gradient), $h_0$는 전엔탈피, $h_w$는 벽면 엔탈피, $\mathrm{Le}$는 루이스 수(Lewis Number), $h_D$는 해리 엔탈피입니다. 정체점 속도 구배는 수정 뉴턴 이론(Modified Newtonian Theory)에 따라 첨두부 곡률 반경 $R_N$에 반비례합니다:

따라서 정체점 열유속은 첨두부 반경의 제곱근에 반비례($q_s \propto 1/\sqrt{R_N}$)합니다. 첨두부가 날카로울수록 항력은 줄어들지만 열유속은 급증하며, 삭마에 의해 첨두부 반경이 커지면($R_N(t) > R_N(0)$) 열유속은 감소하나 형상 항력($C_D$)이 증가하게 됩니다.

2. C/C 복합재 삭마 열화학 반응 및 이동 경계 열전도 지배방정식

탄소-탄소(C/C) 복합재는 $T_w > 1,000 ext{ K}$ 영역에서 확산 제어 산화 반응(Diffusion-controlled Oxidation)을 겪으며, $T_w > 2,800 ext{ K}$ 이상의 초고온에서는 탄소 승화(Sublimation, $C_{(s)} o C_{(g)}, C_2, C_3$)가 지배적으로 일어납니다.

표면 질량 삭마 플럭스 $\dot{m}_{abl}$과 표면 후퇴 속도(Recession Rate) $\dot{s}(t)$는 아레니우스(Arrhenius) 산화 반응 속도론과 산소 확산 한계의 조화 평균으로 모델링됩니다:

삭마로 인해 좌표계의 원점이 후퇴하는 이동 경계(Moving Boundary) 1차원 비정상 열전도 지배방정식은 표면 후퇴 속도 $\dot{s}(t)$에 의한 대류 대류항(Convective-like transport term)을 포함하여 다음과 같이 기술됩니다:

표면 경계 조건($x = 0$)에서의 에너지 보존 법칙은 유입 공력 가열량, 삭마 반응열(흡열 승화 및 발열 산화), 표면 복사 방출열의 평형을 이룹니다:

여기서 $arepsilon$은 C/C 복합재의 방사율(약 0.85~0.92), $\sigma_{SB}$는 슈테판-볼츠만 상수($5.67037 imes 10^{-8} ext{ W/m}^2 ext{K}^4$), $\Delta H_{abl}$은 유효 삭마열(Effective Heat of Ablation, $J/ ext{kg}$)입니다.

3. 첨두부 곡률 후퇴에 따른 공력 계수 변동 및 공력 중심(CP) 드리프트

삭마 과정에서 첨두부 반경 $R_N(t)$가 팽창하면서 탄두의 총 길이 $L(t) = L_0 - s(t)$는 감소하고, 노즈콘 팁 형상이 뭉툭해집니다. 이로 인해 공력 미계수(Aerodynamic Derivatives)가 실시간으로 변동합니다:

공력 중심 $X_{CP}$가 탄두 후방으로 이동하면 정안정성 여유(Static Margin, $ ext{SM} = (X_{CP} - X_{CG})/L_{ref}$)가 과도하게 증가하여 기동성이 저하되고, 반대로 전방으로 이동하면 피칭 불안정성이 급증하여 제어 발산(Loss of Control)을 초래합니다. 따라서 1,000Hz의 온디바이스 유도조종 컴퓨터는 삭마율을 실시간 추정하여 카나드 조종면의 트림 받음각(Trim AoA)과 서보 피드백 게인을 적응형으로 보상해야 합니다.

4. 1,000Hz 0-Alloc C# 극초음속 삭마 열전도 및 공력 보상 비행제어 엔진

아래의 C# 구현체는 탄두 첨두부 삭마 이동 경계 열전도 20개 셀 유한차분법(FDM) 수치해석, Fay-Riddell 열유속, 삭마 후퇴율 $\dot{s}$, 공력 중심 이동량 계산 및 카나드 핀 구동각 명령 산출을 단 하나의 힙 할당도 없이(Zero GC Allocation) 1,000Hz(1.0ms) 주기로 실행하는 고신뢰성 국방 아키텍처입니다.

using System;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using System.Numerics;

namespace ZeliDesk.Defense.Hypersonic
{
    /// <summary>
    /// 탄두 첨두부 형상 및 공력 상태 구조체 (16바이트 정렬)
    /// </summary>
    [StructLayout(LayoutKind.Sequential)]
    public struct NoseGeometryState
    {
        public float NoseRadiusMeters;    // 현재 첨두부 곡률 반경 R_N [m]
        public float TotalRecessionMeters; // 누적 삭마 후퇴량 s(t) [m]
        public float WallTemperatureK;     // 표면 벽면 온도 T_w [K]
        public float HeatFluxWPerM2;       // 표면 열유속 q_s [W/m^2]
        public float DragCoeffCD;          // 삭마 보정 항력계수 C_D
        public float DeltaXcpMeters;       // 공력중심(CP) 이동량 Delta X_CP [m]
        public float PitchTrimDeltaDeg;    // 카나드 피치 트림 보정각 [deg]
    }

    /// <summary>
    /// 극초음속 비행 환경 파라미터 구조체
    /// </summary>
    [StructLayout(LayoutKind.Sequential)]
    public readonly struct HypersonicFlightCondition
    {
        public readonly float MachNumber;       // 마하수 (예: 6.5)
        public readonly float VelocityMps;      // 대기속도 [m/s]
        public readonly float AltitudeMeters;   // 고도 [m]
        public readonly float AirDensityKgM3;   // 대기밀도 rho_inf [kg/m^3]
        public readonly float AmbientTempK;     // 대기온도 T_inf [K]
        public readonly float AngleOfAttackDeg; // 받음각 alpha [deg]

        public HypersonicFlightCondition(float mach, float velocity, float altitude, float density, float ambientTemp, float aoa)
        {
            MachNumber = mach;
            VelocityMps = velocity;
            AltitudeMeters = altitude;
            AirDensityKgM3 = density;
            AmbientTempK = ambientTemp;
            AngleOfAttackDeg = aoa;
        }
    }

    /// <summary>
    /// 1,000Hz 무할당(Zero-Alloc) C/C 복합재 삭마 열전도 및 비행제어 보상 엔진
    /// </summary>
    public ref struct HypersonicAblationEngine
    {
        private const int MeshNodes = 20; // 삭마 열전도 1D 격자점 수
        private const float StefanBoltzmann = 5.67037e-8f; // W/(m^2 K^4)
        private const float Emissivity = 0.88f;            // C/C 복합재 방사율
        private const float DensityCC = 1850.0f;           // C/C 복합재 밀도 [kg/m^3]
        private const float ActivationEnergyR = 21500.0f;  // 활성화 에너지 E_a / R_u [K]
        private const float PreExponentialK0 = 1.25e3f;    // 산화 반응 속도 계수
        private const float HeatOfAblation = 3.2e7f;       // 삭마열 Delta H_abl [J/kg]

        // 1D FDM 온도 분포 버퍼 (스택 할당 Span)
        private Span<float> _temperatureGrid;
        private float _dxMeters;
        private float _initialRadius;
        private float _bodyLength;

        public HypersonicAblationEngine(Span<float> tempBuffer, float totalThickness, float initialRadius, float bodyLength)
        {
            if (tempBuffer.Length < MeshNodes)
                throw new ArgumentException($"Buffer size must be at least {MeshNodes}");

            _temperatureGrid = tempBuffer.Slice(0, MeshNodes);
            _dxMeters = totalThickness / (MeshNodes - 1);
            _initialRadius = initialRadius;
            _bodyLength = bodyLength;
        }

        /// <summary>
        /// 1,000Hz (dt = 0.001s) 주기 비행제어 보상 루프 연산
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public void Update(in HypersonicFlightCondition cond, float dt, ref NoseGeometryState state)
        {
            float currentRn = MathF.Max(0.005f, _initialRadius + state.TotalRecessionMeters);
            state.NoseRadiusMeters = currentRn;

            // 1. Fay-Riddell 정체점 대류 열유속 계산 (qs ~ 1 / sqrt(Rn))
            // qs = C * sqrt(rho_inf / Rn) * V_inf^3
            float rhov3 = MathF.Sqrt(cond.AirDensityKgM3 / currentRn) * MathF.Pow(cond.VelocityMps, 3.0f);
            float qStag = 1.83e-4f * rhov3; // 공칭 열유속 근사치 [W/m^2]
            state.HeatFluxWPerM2 = qStag;

            // 2. 표면 벽면 온도 (격자점 0번)
            float tw = MathF.Max(300.0f, _temperatureGrid[0]);
            state.WallTemperatureK = tw;

            // 3. 아레니우스 산화 반응에 따른 표면 삭마 후퇴율 s_dot [m/s]
            float sDot = 0.0f;
            if (tw > 1100.0f)
            {
                // Arrhenius 삭마 속도론
                float reactionFactor = PreExponentialK0 * MathF.Exp(-ActivationEnergyR / tw);
                sDot = (reactionFactor / DensityCC) * (cond.AirDensityKgM3 * 0.21f);
                sDot = MathF.Min(0.015f, sDot); // 최대 후퇴율 리미터 (15mm/s)
            }

            // 누적 후퇴량 업데이트
            state.TotalRecessionMeters += sDot * dt;

            // 4. 이동 경계 1D 열전도 FDM 수치적분
            SolveThermalGrid(qStag, tw, sDot, dt);

            // 5. 첨두부 반경 팽창에 따른 공력계수 변화량 산출
            // 형상 뭉툭화 -> 항력 증가, 양력 감소, 공력중심 후방 이동
            float bluntRatio = state.TotalRecessionMeters / _initialRadius;
            state.DragCoeffCD = 0.18f + 0.12f * MathF.Tanh(bluntRatio * 2.0f);

            // 공력중심 이동량 Delta X_CP: 양의 방향은 후방 이동 [m]
            float deltaXcp = 0.085f * MathF.Tanh(bluntRatio * 1.5f);
            state.DeltaXcpMeters = deltaXcp;

            // 6. 공력중심 이동에 따른 카나드 피치 트림 보정 제어각 산출
            // 모멘트 평형 조건: Delta Cm_abl = Delta X_cp * C_L / c_bar
            float dynamicPressure = 0.5f * cond.AirDensityKgM3 * cond.VelocityMps * cond.VelocityMps;
            float aeroFeedbackGain = 1.45f * (dynamicPressure / 1.5e6f);
            float trimCorrectionDeg = -(deltaXcp / _bodyLength) * 57.2958f * aeroFeedbackGain;

            state.PitchTrimDeltaDeg = Math.Clamp(trimCorrectionDeg, -8.0f, 8.0f);
        }

        /// <summary>
        /// 이동 경계 열전도 FDM 유한차분 해석 (Zero-Alloc 스택 연산)
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private void SolveThermalGrid(float qStag, float tw, float sDot, float dt)
        {
            // C/C 복합재 열물성 (온도 종속 근사)
            float kTh = 35.0f; // 열전도율 [W/(m K)]
            float cp = 1400.0f; // 비열 [J/(kg K)]
            float alpha = kTh / (DensityCC * cp); // 열확산율 [m^2/s]

            // 복사 방출열 손실 및 삭마 흡열량
            float qRad = Emissivity * StefanBoltzmann * MathF.Pow(tw, 4.0f);
            float qAbl = DensityCC * sDot * HeatOfAblation;
            float qNetSurface = qStag - qRad - qAbl;

            // 경계조건 0번 표면 노드 온도 업데이트 (전진 오일러)
            float dT0_dx = ( _temperatureGrid[1] - _temperatureGrid[0] ) / _dxMeters;
            float surfaceFluxTerm = (qNetSurface + kTh * dT0_dx) / (DensityCC * cp * (_dxMeters * 0.5f));
            _temperatureGrid[0] += surfaceFluxTerm * dt;

            // 내부 노드 FDM (대류 대류항 sDot * dT/dx 반영)
            Span<float> nextTemp = stackalloc float[MeshNodes];
            nextTemp[0] = _temperatureGrid[0];

            for (int i = 1; i < MeshNodes - 1; i++)
            {
                float d2T_dx2 = (_temperatureGrid[i + 1] - 2.0f * _temperatureGrid[i] + _temperatureGrid[i - 1]) / (_dxMeters * _dxMeters);
                float dT_dx = (_temperatureGrid[i + 1] - _temperatureGrid[i - 1]) / (2.0f * _dxMeters);

                // 이동 경계 대류항: sDot * dT/dx
                float dT_dt = alpha * d2T_dx2 + sDot * dT_dx;
                nextTemp[i] = _temperatureGrid[i] + dT_dt * dt;
            }

            // 후방 경계 절열 조건 (dT/dx = 0)
            nextTemp[MeshNodes - 1] = nextTemp[MeshNodes - 2];

            // 버퍼 갱신
            for (int i = 0; i < MeshNodes; i++)
            {
                _temperatureGrid[i] = nextTemp[i];
            }
        }
    }
}

5. 2022 개정 공학·물리 교육과정 및 국방 실무 연계