초음속 제트엔진 터빈 블레이드 단결정 초합금 고온 크리프(Creep) & 1,000Hz 0-Alloc C# 라르손-밀러 수명 예측 엔진

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Stem.Propulsion;

public readonly ref struct TurbineBladeHealthState
{
    public readonly double CentrifugalStressMpa;
    public readonly double CreepStrainRatePerHour;
    public readonly double AccumulatedDamageFraction;
    public readonly double RemainingLifeHours;
    public readonly bool IsReplacementRequired;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public TurbineBladeHealthState(double stress, double strainRate, double damage, double remLife, bool replace)
    {
        CentrifugalStressMpa = stress;
        CreepStrainRatePerHour = strainRate;
        AccumulatedDamageFraction = damage;
        RemainingLifeHours = remLife;
        IsReplacementRequired = replace;
    }
}

public ref struct TurbineBladeCreepSimulator
{
    private const double DensityKgM3 = 8700.0; // CMSX-4 밀도
    private const double BladeLengthM = 0.08; // 80 mm 블레이드
    private const double RotorRadiusM = 0.35; // 350 mm 로터 반경
    private const double ConstantC = 20.0; // 라르손-밀러 상수

    private double _accumulatedDamage;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public TurbineBladeCreepSimulator(double initialDamage = 0.05)
    {
        _accumulatedDamage = initialDamage;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public TurbineBladeHealthState Step(double rotorRpm, double gasTempCelsius, double dtSeconds = 0.001)
    {
        // 1. 블레이드 루트 원심 인장 응력 산출 (MPa)
        double omegaRadS = rotorRpm * (2.0 * Math.PI / 60.0);
        double rMean = RotorRadiusM + BladeLengthM * 0.5;
        double stressPa = DensityKgM3 * omegaRadS * omegaRadS * rMean * BladeLengthM * 0.65; // 테이퍼 단면 계수
        double stressMpa = stressPa * 1e-6;

        // 2. 블레이드 금속 유효 온도 산출 (내부 냉각 홀 효과 고려: 가스온도 대비 ~300K 저감)
        double metalTempK = (gasTempCelsius + 273.15) - 280.0;

        // 3. 라르손-밀러 파라미터 (LMP) 계산: LMP = a0 - a1 * log10(sigma)
        double logStress = Math.Log10(Math.Max(10.0, stressMpa));
        double lmp = 32.5 - 5.2 * logStress; // 10^3 단위

        // 파단 수명 시간 (hours)
        double exponent = (lmp * 1000.0 / metalTempK) - ConstantC;
        double ruptureLifeHours = Math.Pow(10.0, Math.Min(8.0, Math.Max(-1.0, exponent)));

        // 4. 노턴 정상 크리프 변형률 속도 (%/hour)
        double strainRatePerHour = 1.2e-5 * Math.Pow(stressMpa / 150.0, 5.5) * Math.Exp(-32000.0 / metalTempK);

        // 5. 마이너 누적 손상 적분 (dt = 0.001초 -> hours 변환)
        double dtHours = dtSeconds / 3600.0;
        double incrementalDamage = dtHours / ruptureLifeHours;
        _accumulatedDamage += incrementalDamage;

        double remainingLifeHours = Math.Max(0.0, (1.0 - _accumulatedDamage) * ruptureLifeHours);
        bool replaceRequired = _accumulatedDamage >= 0.85;

        return new TurbineBladeHealthState(stressMpa, strainRatePerHour, _accumulatedDamage, remainingLifeHours, replaceRequired);
    }
}