초소형 가스터빈(MGT) 공기 포일 베어링(Air Foil Bearing) 고속 로터 동역학 & 1,000Hz 0-Alloc C# 유체동압 틸팅패드 서보 엔진

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Stem.Propulsion.AirFoilBearing
{
    /// <summary>
    /// 공기 포일 베어링 센서 계측 및 운전 상태 구조체 (0-Alloc Stack Only)
    /// </summary>
    public readonly ref struct BearingTelemetry
    {
        public readonly float RotorRpm;          // [RPM] 로터 회전수 (예: 120,000)
        public readonly float DisplacementXUm;   // [um] X축 축 변위
        public readonly float DisplacementYUm;   // [um] Y축 축 변위
        public readonly float BearingTempC;      // [deg C] 베어링 온도
        public readonly float AmbientPressurePa; // [Pa] 흡입 대기압 (예: 101325)
        public readonly float BearingRadiusM;    // [m] 저널 반경 (예: 0.015m)
        public readonly float BearingLengthM;    // [m] 베어링 축 길이 (예: 0.030m)
        public readonly float RadialClearanceUm; // [um] 공칭 반경 틈새 (예: 25um)
        public readonly float BumpStiffnessNm3;  // [N/m^3] 범프 포일 단위면적 구조강성

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public BearingTelemetry(
            float rotorRpm,
            float displacementXUm,
            float displacementYUm,
            float bearingTempC,
            float ambientPressurePa,
            float bearingRadiusM,
            float bearingLengthM,
            float radialClearanceUm,
            float bumpStiffnessNm3)
        {
            RotorRpm = rotorRpm;
            DisplacementXUm = displacementXUm;
            DisplacementYUm = displacementYUm;
            BearingTempC = bearingTempC;
            AmbientPressurePa = ambientPressurePa;
            BearingRadiusM = bearingRadiusM;
            BearingLengthM = bearingLengthM;
            RadialClearanceUm = radialClearanceUm;
            BumpStiffnessNm3 = bumpStiffnessNm3;
        }
    }

    /// <summary>
    /// 공기 포일 베어링 서보 연산 결과 구조체
    /// </summary>
    public readonly ref struct BearingAnalysisResult
    {
        public readonly float DynamicViscosityPaS;  // [Pa*s] 고온 공기 점성계수
        public readonly float BearingNumberLambda;  // [무차원] 압축성 베어링 수
        public readonly float EccentricityRatio;    // [0~1] 로터 편심률 (epsilon)
        public readonly float EstimatedLoadN;       // [N] 유체동압 지지 하중 용량
        public readonly float MinimumFilmThicknessUm; // [um] 최소 유막 두께 (h_min)
        public readonly bool IsLiftedOff;           // 기체막 완전 부양 여부
        public readonly bool IsWhirlUnstable;       // 서브동기 진동 불안정 경보

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public BearingAnalysisResult(
            float dynamicViscosityPaS,
            float bearingNumberLambda,
            float eccentricityRatio,
            float estimatedLoadN,
            float minimumFilmThicknessUm,
            bool isLiftedOff,
            bool isWhirlUnstable)
        {
            DynamicViscosityPaS = dynamicViscosityPaS;
            BearingNumberLambda = bearingNumberLambda;
            EccentricityRatio = eccentricityRatio;
            EstimatedLoadN = estimatedLoadN;
            MinimumFilmThicknessUm = minimumFilmThicknessUm;
            IsLiftedOff = isLiftedOff;
            IsWhirlUnstable = isWhirlUnstable;
        }
    }

    /// <summary>
    /// 1,000Hz 0-Alloc 초소형 가스터빈 공기 포일 베어링 서보 엔진
    /// </summary>
    public static class AirFoilBearingServoEngine
    {
        private const float MinSafeFilmThicknessUm = 3.5f; // [um] 최소 안전 유막 두께
        private const float MaxStableEccentricity = 0.85f;  // 최대 허용 편심률

        /// <summary>
        /// 1ms 주기 실시간 유체동압 유막 형성 및 진동 안정성 해석 수행
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public static BearingAnalysisResult EvaluateBearing(in BearingTelemetry telemetry)
        {
            // 1. 고온 공기 점성계수 계산 (Sutherland 공식)
            // mu(T) = mu_0 * (T/T_0)^(3/2) * (T_0 + S) / (T + S)
            float tempKelvin = telemetry.BearingTempC + 273.15f;
            const float t0 = 273.15f;
            const float mu0 = 1.716e-5f; // [Pa*s]
            const float sConst = 110.4f; // [K]
            float sutherlandRatio = (tempKelvin / t0);
            float visFactor = (t0 + sConst) / (tempKelvin + sConst);
            float muAir = mu0 * MathF.Pow(sutherlandRatio, 1.5f) * visFactor;

            // 2. 로터 각속도 및 무차원 베어링 수(Compressibility Bearing Number Lambda)
            float omega = telemetry.RotorRpm * (2.0f * MathF.PI / 60.0f);
            float clearanceMeters = telemetry.RadialClearanceUm * 1e-6f;
            float rJournal = telemetry.BearingRadiusM;
            float pAtm = telemetry.AmbientPressurePa;

            float lambda = (6.0f * muAir * omega * rJournal * rJournal) / 
                           (pAtm * clearanceMeters * clearanceMeters);

            // 3. 로터 편심률(Eccentricity Ratio) 계산
            float orbitRadiusUm = MathF.Sqrt(
                telemetry.DisplacementXUm * telemetry.DisplacementXUm + 
                telemetry.DisplacementYUm * telemetry.DisplacementYUm);
            float epsilon = orbitRadiusUm / telemetry.RadialClearanceUm;
            if (epsilon > 0.98f) epsilon = 0.98f;

            // 4. 범프 포일 탄성 변형을 반영한 유효 최소 유막 두께 계산
            // h_min = c * (1 - epsilon) + delta_bump
            float rigidHMinUm = telemetry.RadialClearanceUm * (1.0f - epsilon);
            // 압력 상승에 의한 범프 변형 근사 (포일 처짐으로 인한 틈새 확장 보정)
            float bumpDeflectionUm = (lambda * 0.15f) * (epsilon / (1.0f - epsilon * epsilon + 0.05f));
            float effectiveHMinUm = rigidHMinUm + bumpDeflectionUm;

            // 5. 유체동압 지지 부하 용량(Load Capacity W) 산출
            // W = p_a * R * L * (pi/2) * lambda * epsilon / sqrt(1 - epsilon^2) (단축 베어링 근사치)
            float geoFactor = pAtm * rJournal * telemetry.BearingLengthM;
            float denom = MathF.Sqrt(MathF.Max(0.01f, 1.0f - epsilon * epsilon));
            float loadCapacityN = geoFactor * (MathF.PI * 0.5f) * (lambda / (1.0f + 0.5f * lambda)) * (epsilon / denom);

            // 6. 이탈(Lift-off) 상태 판별: 최소 유막 두께가 거칠기 한계 이상이고 회전수 기준 충족 시
            bool isLifted = (effectiveHMinUm > MinSafeFilmThicknessUm) && (telemetry.RotorRpm > 15000.0f);

            // 7. 서브동기 진동 훨(Whirl) 불안정성 판별: 편심률 과다 또는 임계 베어링 수 초과 시
            bool isWhirlUnstable = (epsilon > MaxStableEccentricity) || 
                                   (orbitRadiusUm > telemetry.RadialClearanceUm * 0.7f && telemetry.RotorRpm > 80000.0f);

            return new BearingAnalysisResult(
                dynamicViscosityPaS: muAir,
                bearingNumberLambda: lambda,
                eccentricityRatio: epsilon,
                estimatedLoadN: loadCapacityN,
                minimumFilmThicknessUm: effectiveHMinUm,
                isLiftedOff: isLifted,
                isWhirlUnstable: isWhirlUnstable
            );
        }
    }
}