대한민국 공군 중고도 정찰무인기(KUS-FS) 동계 20,000ft 결빙(Icing) 조우 시 전열 방빙 패드 듀티비 제어 및 1,000Hz 0-Alloc C# 비행 안정성 복원 엔진

1. 과냉각 수적(Supercooled Droplets) 충돌 역학 및 주익 포집 계수($\beta$)

대기 중 액체 수분 함량(Liquid Water Content, $\text{LWC}, \text{g/m}^3$)과 수적 중위 체적 직경(Median Volume Diameter, $\text{MVD}, \mu\text{m}$)에 따른 국소 물방울 충돌 질량 유속 $\dot{m}_{imp}(s)$는 다음과 같이 기술됩니다:

여기서 $s$는 주익 전연 표면 좌표, $\beta(s) = dy_0 / ds$는 공기역학적 유선(Streamline)을 이탈하여 에어포일에 직접 충돌하는 수적 궤적의 국소 포집 효율(Local Collection Efficiency)입니다. 수정 관성 파라미터 $K_0$와 레이놀즈수 $\text{Re}_d$에 따라 수적의 관성이 클수록($\text{MVD} > 20\mu\text{m}$) 유선을 벗어나 전연 정체점 부근에 집중 충돌합니다.

2. Messinger 열평형 방정식 및 전열 방빙 요구 열유속

결빙을 방지하고 충돌 수적을 즉각 기화 또는 미결빙 유출(Runback Water) 상태로 유지하기 위한 주익 표면 미소 면적에서의 Messinger 에너지 보존 방정식은 다음과 같습니다:

여기서 $q_{in}$은 전열 방빙 패드에서 공급되는 단위면적당 열유속($\text{W/m}^2$)이며, 주요 항들은 다음과 같이 계산됩니다:

완전 증발 방빙(Evaporative Anti-icing) 조건에서는 $q_{in} \ge q_{conv} + q_{evap} + q_{sens} - q_{aero} - q_{kinetic}$을 만족해야 하므로, 외기온이 $-20^\circ\text{C}$이고 풍속이 $80\text{ m/s}$일 때 약 $15 \sim 25\text{ kW/m}^2$의 전력이 요구됩니다.

3. 결빙층 성장에 따른 양력 계수 저하 및 공력 중심(CP) 변동

방빙 전력이 부족하여 리딩에지에 미세 결빙층(두께 $t_{ice}$)이 형성될 경우, 흡입면 압력 피크가 소멸하고 조기 경계층 박리가 발생합니다:

공력 중심 $X_{CP}$가 전방으로 이동하면서 피칭 업 모멘트가 유발되고 실속 받음각($\alpha_{stall}$)이 기존 $14^\circ$에서 $8^\circ$ 수준으로 급격히 축소됩니다. 따라서 비행제어 컴퓨터(FCC)는 결빙 감지 즉시 엘리베이터 하향 트림각을 인가하고 자동 엔진 스로틀을 보상해야 합니다.

4. 1,000Hz 0-Alloc C# 전열 방빙 듀티비 제어 및 비행안정성 복원 엔진

아래의 C# 구현체는 외기온, 대기속도, LWC, 광학 결빙 센서 신호를 입력받아 주익 4개 존의 PWM 전열 방빙 듀티비와 결빙 양력 손실을 보상하는 엘리베이터 트림각을 가비지 컬렉션 부하 0%로 1,000Hz 주기로 연산하는 고신뢰성 무인기 비행제어 코어입니다.

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

namespace ZeliDesk.Defense.UavAvionics
{
    /// <summary>
    /// 비행 환경 및 결빙 기상 관측 데이터
    /// </summary>
    [StructLayout(LayoutKind.Sequential)]
    public readonly struct IcingAtmosphereFrame
    {
        public readonly float TrueAirspeedMps;   // 진대기속도 TAS [m/s] (예: 75.0)
        public readonly float AmbientTempC;      // 외기온도 T_inf [C] (예: -18.0)
        public readonly float AmbientPressurePa; // 정압력 P_inf [Pa] (고도 20,000ft: 약 46,500)
        public readonly float LiquidWaterContent;// LWC [g/m^3] (예: 0.45)
        public readonly float IceThicknessMm;    // 측정된 주익 결빙 두께 [mm]
        public readonly float AngleOfAttackDeg;  // 현재 받음각 [deg]

        public IcingAtmosphereFrame(float tas, float temp, float press, float lwc, float iceThick, float aoa)
        {
            TrueAirspeedMps = tas;
            AmbientTempC = temp;
            AmbientPressurePa = press;
            LiquidWaterContent = lwc;
            IceThicknessMm = iceThick;
            AngleOfAttackDeg = aoa;
        }
    }

    /// <summary>
    /// 전열 방빙 히터 PWM 지령 및 비행제어 보상 지령
    /// </summary>
    [StructLayout(LayoutKind.Sequential)]
    public struct AntiIcingControlCommand
    {
        public float LeadingEdgeDutyCycle; // 전연 정체점 존 PWM 듀티비 (0.0~1.0)
        public float UpperSurgDutyCycle;   // 상면 흡입면 존 PWM 듀티비 (0.0~1.0)
        public float LowerSurgDutyCycle;   // 하면 압력면 존 PWM 듀티비 (0.0~1.0)
        public float TotalHeaterPowerKw;   // 소모 총 방빙 전력 [kW]
        public float ElevatorTrimOffsetDeg;// 결빙 모멘트 보상 엘리베이터 트림각 [deg]
        public float ThrottleBonusPercent; // 항력 증가 보상 추가 엔진 스로틀 [%]
        public bool IcingStallWarning;     // 실속 여유 위험 경고
    }

    /// <summary>
    /// 1,000Hz 무할당(Zero-Alloc) KUS-FS 전열 방빙 및 비행안정성 복원 엔진
    /// </summary>
    public ref struct KusFsAntiIcingFlightEngine
    {
        private const float WingAreaM2 = 26.5f;          // KUS-FS 주익 면적
        private const float MeanAerodynamicChord = 1.45f; // 평균 공력 시위선 c [m]
        private const float MaxHeaterPowerKw = 18.0f;    // 발전기 가용 최대 방빙 전력 (18 kW)
        private const float AirCp = 1005.0f;             // 공기 정압비열 [J/(kg K)]
        private const float LatentHeatVapor = 2.5e6f;    // 기화잠열 [J/kg]

        public KusFsAntiIcingFlightEngine()
        {
        }

        /// <summary>
        /// 1,000Hz (dt = 0.001s) 주기 방빙 전력 최적화 및 비행제어 보상 루프
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public void Update(in IcingAtmosphereFrame env, ref AntiIcingControlCommand cmd)
        {
            // 1. 대류 열전달 계수 h_c 추정 (평판 난류 열전달 근사 h_c ~ V^0.8)
            float airDensity = env.AmbientPressurePa / (287.05f * (env.AmbientTempC + 273.15f));
            float dynViscosity = 1.65e-5f;
            float reynolds = (airDensity * env.TrueAirspeedMps * MeanAerodynamicChord) / dynViscosity;
            float nusselt = 0.037f * MathF.Pow(reynolds, 0.8f) * MathF.Pow(0.71f, 0.33f);
            float thermConductivity = 0.024f;
            float hc = (nusselt * thermConductivity) / MeanAerodynamicChord;

            // 2. Messinger 열평형 기반 단위면적당 요구 열유속 계산 [W/m^2]
            // 목표 표면온도 T_w = +3.0 C (미결빙 방빙)
            float targetTw = 3.0f;
            float deltaT = targetTw - env.AmbientTempC;

            float qConv = hc * deltaT;
            float mImp = env.TrueAirspeedMps * (env.LiquidWaterContent * 0.001f) * 0.65f; // kg/(m^2 s)
            float qSens = mImp * 4186.0f * deltaT;
            float qEvap = (0.622f * hc / (AirCp * env.AmbientPressurePa)) * LatentHeatVapor * 450.0f;
            float qAeroKinetic = (0.85f * (env.TrueAirspeedMps * env.TrueAirspeedMps) / (2.0f * AirCp) * hc) + (0.5f * mImp * env.TrueAirspeedMps * env.TrueAirspeedMps);

            float qNetRequired = MathF.Max(0.0f, qConv + qSens + qEvap - qAeroKinetic); // [W/m^2]

            // 3. 발전기 가용 전력 한도 내 전열 매트 듀티비 최적 분배
            // 전연(Leading Edge) 60%, 상면(Upper) 25%, 하면(Lower) 15%
            float totalRequiredKw = (qNetRequired * WingAreaM2 * 0.12f) / 1000.0f; // 방빙 면적 약 12%
            float powerScale = (totalRequiredKw > MaxHeaterPowerKw) ? (MaxHeaterPowerKw / totalRequiredKw) : 1.0f;
            cmd.TotalHeaterPowerKw = MathF.Min(totalRequiredKw, MaxHeaterPowerKw);

            cmd.LeadingEdgeDutyCycle = Math.Clamp((totalRequiredKw / MaxHeaterPowerKw) * powerScale * 1.0f, 0.0f, 1.0f);
            cmd.UpperSurgDutyCycle = Math.Clamp(cmd.LeadingEdgeDutyCycle * 0.75f, 0.0f, 1.0f);
            cmd.LowerSurgDutyCycle = Math.Clamp(cmd.LeadingEdgeDutyCycle * 0.50f, 0.0f, 1.0f);

            // 4. 결빙층 형성에 따른 공력 손실 및 실속 위험 평가
            float iceRatio = env.IceThicknessMm / (MeanAerodynamicChord * 1000.0f);
            float stallAoADeg = 14.5f - (35.0f * MathF.Pow(iceRatio, 0.5f));
            cmd.IcingStallWarning = (env.AngleOfAttackDeg > (stallAoADeg - 2.5f));

            // 5. 엘리베이터 트림각 및 스로틀 보상 지령 산출
            // 결빙 시 양력 감소 및 항력 증가 -> 기수 내림 트림 및 엔진 스로틀 보상
            float deltaElevator = -1.8f * MathF.Sqrt(env.IceThicknessMm); // 기수 유지용 하향 편향각
            cmd.ElevatorTrimOffsetDeg = Math.Clamp(deltaElevator, -6.0f, 0.0f);

            float deltaThrottle = 3.5f * env.IceThicknessMm; // 항력 증가 극복용 스로틀 [%]
            cmd.ThrottleBonusPercent = Math.Clamp(deltaThrottle, 0.0f, 25.0f);
        }
    }
}

5. 2022 개정 물리·지구과학·정보 교육과정 연계표