초전도 양자컴퓨팅 극저온 희석냉동기(Dilution Refrigerator) 3He-4He 상분리 & 1,000Hz 0-Alloc C# 10mK 열평형 제어 엔진

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Stem.Cryogenics;

public readonly ref struct DilutionFridgeStatus
{
    public readonly double MixingChamberTempMilliKelvin;
    public readonly double CoolingPowerMicroWatts;
    public readonly double CirculationRateMicromolSec;
    public readonly double HeaterPowerMicroWatts;
    public readonly bool IsQpuReady;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public DilutionFridgeStatus(double tempMk, double coolUv, double flowUmol, double heaterUv, bool ready)
    {
        MixingChamberTempMilliKelvin = tempMk;
        CoolingPowerMicroWatts = coolUv;
        CirculationRateMicromolSec = flowUmol;
        HeaterPowerMicroWatts = heaterUv;
        IsQpuReady = ready;
    }
}

public ref struct DilutionRefrigeratorSimulator
{
    private double _mixingChamberTempK;
    private double _stillTempK;
    private double _heaterPowerWatts;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public DilutionRefrigeratorSimulator(double initialTempK = 0.050) // 초기 50 mK
    {
        _mixingChamberTempK = initialTempK;
        _stillTempK = 0.85; // Still 정상 운전 0.85 K
        _heaterPowerWatts = 1.0e-3; // 1 mW
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public DilutionFridgeStatus Step(double qpuHeatLoadWatts, double dtSeconds = 0.001)
    {
        // 1. Still 증기압 기반 3He 순환율 계산 (Antoine 증기압 근사)
        // 0.85K에서 3He 증기압이 4He보다 수백 배 높음
        double circulationMolPerSec = 2.0e-4 * Math.Exp(-0.5 / _stillTempK);

        // 2. 혼합실 유효 냉각능 산출 (J/s = W)
        // Tin은 열교환기를 통과하여 대략 Tmc의 1.3배 수준으로 예냉되었다고 가정
        double tIn = _mixingChamberTempK * 1.3;
        double qCoolWatts = circulationMolPerSec * (84.0 * _mixingChamberTempK * _mixingChamberTempK - 12.0 * tIn * tIn);
        if (qCoolWatts < 0.0) qCoolWatts = 0.0;

        // 3. 기생 열부하 + QPU 마이크로파 발열
        double parasiticLoadWatts = 0.5e-6; // 0.5 uW 차폐 잔류열
        double netHeatWatts = parasiticLoadWatts + qpuHeatLoadWatts - qCoolWatts;

        // 4. 열용량(C_v = alpha * T) 및 온도 미분 방정식 전진 적분
        double thermalMass = 0.02 * _mixingChamberTempK; // J/K
        double dT = (netHeatWatts / Math.Max(1e-6, thermalMass)) * dtSeconds;
        _mixingChamberTempK += dT;
        if (_mixingChamberTempK < 0.005) _mixingChamberTempK = 0.005; // 5 mK 물리적 하한

        // 5. 목표 10 mK 유지를 위한 Still 히터 능동 PID 제어
        double targetK = 0.010;
        double error = _mixingChamberTempK - targetK;
        if (error > 0.0)
        {
            _stillTempK += error * 0.05 * dtSeconds; // Still 온도 올려 순환율 증강
        }
        else
        {
            _stillTempK -= Math.Abs(error) * 0.05 * dtSeconds;
        }
        if (_stillTempK < 0.7) _stillTempK = 0.7;
        if (_stillTempK > 1.0) _stillTempK = 1.0;

        double tempMk = _mixingChamberTempK * 1000.0;
        double coolUw = qCoolWatts * 1.0e6;
        double flowUmol = circulationMolPerSec * 1.0e6;
        double heaterUw = _heaterPowerWatts * 1.0e6;
        bool ready = tempMk <= 15.0; // 15 mK 이하 시 양자 연산 준비 완료

        return new DilutionFridgeStatus(tempMk, coolUw, flowUmol, heaterUw, ready);
    }
}