극초음속 비행체 능동 열방호 발한 냉각(Transpiration Cooling) & 1,000Hz 0-Alloc C# 다공성 매질 이상 유동 서보 엔진

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Stem.HypersonicThermal
{
    /// <summary>
    /// 극초음속 비행체 다공성 매질 발한 냉각 이상유동 서보 상태
    /// 0-Alloc ref struct: GC 압압 완전 배제 (1,000Hz 실시간 루프)
    /// </summary>
    public ref struct TranspirationCoolingState
    {
        public float WallTemperatureK;       // 표면 벽면 온도 [K]
        public float TargetTemperatureK;     // 목표 안전 온도 [K] (예: 1100 K)
        public float AeroHeatFluxWm2;        // 외부 극초음속 공력 열유속 [W/m^2]
        public float BoilingFrontDepthM;     // 기화 전방 깊이 [m] (0 = 표면, 음수 = 내부)
        public float CoolantPressurePa;      // 냉각재 공급 압력 [Pa]
        public float MassFluxKgM2s;          // 단위면적당 질량 유속 [kg/(m^2*s)]
        public float PumpDutyCycle;          // 냉각재 고압 펌프 PWM 제어 듀티비 (0.0 ~ 1.0)
        public bool IsDryoutWarning;         // 건조 파괴 위험 플래그
    }

    /// <summary>
    /// 1,000Hz 무할당 다공성 발한냉각 이상유동 적응형 서보 엔진
    /// </summary>
    public static class TranspirationCoolingServoEngine
    {
        // 다공성 흑연/세라믹 물성 상수
        private const float Porosity = 0.35f;               // 기공률 epsilon
        private const float PoreRadiusM = 5.0e-6f;          // 유효 기공 반경 [m]
        private const float PermeabilityM2 = 2.5e-12f;      // 고유 침투율 K [m^2]
        private const float ForchheimerBeta = 1.2e6f;       // 비선형 관성 계수 beta [1/m]
        private const float WallThicknessM = 0.015f;        // 다공체 두께 [m] (15mm)

        // 냉각재(초순수) 열역학 상수
        private const float LatentHeatHfg = 2.257e6f;       // 기화 잠열 [J/kg]
        private const float CpLiquid = 4184.0f;             // 액체 비열 [J/(kg*K)]
        private const float CpVapor = 2010.0f;              // 증기 비열 [J/(kg*K)]
        private const float SaturationTempK = 373.15f;      // 포화 온도 (1 atm 기준)
        private const float LiquidDensity = 960.0f;         // 평균 액체 밀도 [kg/m^3]
        private const float LiquidViscosity = 2.8e-4f;      // 액체 점도 [Pa*s]

        /// <summary>
        /// 1,000Hz (1ms) 단일 실행 프레임: 공력 열유속 및 벽면 온도 기반 냉각재 유량 갱신
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static void UpdateFrame(
            ref TranspirationCoolingState state,
            float measuredWallTempK,
            float currentHeatFluxWm2,
            float dtSeconds)
        {
            state.WallTemperatureK = measuredWallTempK;
            state.AeroHeatFluxWm2 = currentHeatFluxWm2;

            // 1) 열유속 방어에 요구되는 이론 최소 질량 유속 계산 (잠열 + 현열 + 불어내기 효과)
            float deltaT_liquid = MathF.Max(0.0f, SaturationTempK - 300.0f);
            float deltaT_vapor = MathF.Max(0.0f, state.WallTemperatureK - SaturationTempK);
            float effectiveEnthalpy = LatentHeatHfg + (CpLiquid * deltaT_liquid) + (CpVapor * deltaT_vapor);

            // 불어내기 효과(Blowing Effect) 근사: 고온 벽면 부근 경계층 박리 차폐 계수
            float h0 = 1500.0f; // 기본 대류 열전달 계수 [W/(m^2*K)]
            float bFactor = (state.MassFluxKgM2s * CpVapor) / (h0 + 1e-4f);
            float blowingReduction = bFactor > 0.01f ? MathF.Log(1.0f + 2.0f * bFactor) / (2.0f * bFactor) : 1.0f;
            float netAeroFlux = currentHeatFluxWm2 * blowingReduction;

            float targetMassFlux = netAeroFlux / (effectiveEnthalpy + 1e-3f);

            // 2) 벽면 온도 오차 기반 PID 피드백 결합 (목표 온도 유지)
            float tempError = state.WallTemperatureK - state.TargetTemperatureK;
            float kp = 0.0025f;
            float feedbackFlux = kp * tempError;
            float requiredFlux = MathF.Max(0.05f, targetMassFlux + feedbackFlux);

            // 3) Ergun-Forchheimer 다공성 압력 손실 모델 역산 -> 요구 공급 압력 산출
            float u = requiredFlux / LiquidDensity; // 여과 속도
            float dP_viscous = (LiquidViscosity / PermeabilityM2) * u * WallThicknessM;
            float dP_inertial = ForchheimerBeta * LiquidDensity * (u * u) * WallThicknessM;
            float requiredPressure = 101325.0f + dP_viscous + dP_inertial; // 절대 압력 [Pa]

            // 4) 모세관 한계 및 건조(Dryout) 여부 판정
            float maxCapillaryPressure = (2.0f * 0.058f) / PoreRadiusM; // sigma = 0.058 N/m
            state.IsDryoutWarning = (dP_viscous + dP_inertial) > (maxCapillaryPressure * 1.5f) && (state.WallTemperatureK > 1400.0f);

            // 5) 펌프 듀티비 변조 (0.0 ~ 1.0)
            float maxPumpPressure = 2.5e6f; // 최대 2.5 MPa (25 bar)
            float duty = MathF.Min(1.0f, MathF.Max(0.05f, requiredPressure / maxPumpPressure));

            state.MassFluxKgM2s = requiredFlux;
            state.CoolantPressurePa = requiredPressure;
            state.PumpDutyCycle = duty;

            // 기화 전방(Boiling Front) 위치 추정: 표면(0)으로부터 내부 깊이 [m]
            if (state.WallTemperatureK > SaturationTempK)
            {
                float excessTemp = state.WallTemperatureK - SaturationTempK;
                state.BoilingFrontDepthM = MathF.Min(0.0f, -0.00001f * excessTemp);
            }
            else
            {
                state.BoilingFrontDepthM = 0.0f;
            }
        }
    }
}