초전도 토카막 핵융합로 디버터(Divertor) 텅스텐 플라즈마 대면재(PFC) 초고열유속(10MW/m²) 열응력 수치해석 & 1,000Hz 0-Alloc C# 유한차분 열전도 분산 냉각 제어 엔진

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Stem.Fusion;

public readonly ref struct DivertorStatus
{
    public readonly double SurfaceTemperature;
    public readonly double CoolingInterfaceTemperature;
    public readonly double MaxThermalStressMpa;
    public readonly double RequiredCoolantFlowKgS;
    public readonly bool IsRecrystallizationRisk;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public DivertorStatus(double surfaceTemp, double coolingTemp, double stressMpa, double flowKgS, bool risk)
    {
        SurfaceTemperature = surfaceTemp;
        CoolingInterfaceTemperature = coolingTemp;
        MaxThermalStressMpa = stressMpa;
        RequiredCoolantFlowKgS = flowKgS;
        IsRecrystallizationRisk = risk;
    }
}

public ref struct DivertorThermalSimulator
{
    private const int GridNodes = 64; // 공간 이산화 격자수
    private const double BlockThickness = 0.025; // 25mm 텅스텐 + 구리 튜브
    private const double Dx = BlockThickness / (GridNodes - 1);
    private const double Dt = 0.001; // 1ms 타임스텝 (1,000Hz)

    // 재료 물성치 (텅스텐 평균)
    private const double DensityW = 19250.0; // kg/m^3
    private const double SpecificHeatW = 138.0; // J/(kg*K)
    private const double ThermalCondW = 160.0; // W/(m*K)
    private const double ElasticModulus = 400.0e9; // 400 GPa
    private const double ThermalExpansion = 4.5e-6; // 1/K
    private const double PoissonRatio = 0.28;
    private const double RecrystallizationLimit = 1300.0; // °C

    // 스택 격자 온도 배열 (0-Alloc)
    private Span<double> _temperatures;
    private double _coolantTemp;
    private double _currentCoolantFlow;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public DivertorThermalSimulator(Span<double> nodeBuffer, double initialTemp = 150.0)
    {
        if (nodeBuffer.Length < GridNodes)
            throw new ArgumentException("Buffer must have at least 64 nodes.");

        _temperatures = nodeBuffer.Slice(0, GridNodes);
        _temperatures.Fill(initialTemp);
        _coolantTemp = 120.0; // 가압수 120°C
        _currentCoolantFlow = 5.0; // 기본 5 kg/s
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public DivertorStatus Step(double plasmaHeatFluxWattsM2)
    {
        // 1. 대류 열전달 계수 산출 (Dittus-Boelter 난류 상관식 기반)
        double hCool = 25000.0 + (_currentCoolantFlow * 4000.0); // W/(m^2*K)

        // 2. 유한차분법(FDM) 전진 오일러 / 확산 스킴
        double alpha = ThermalCondW / (DensityW * SpecificHeatW);
        double r = alpha * Dt / (Dx * Dx);

        // 표면 노드 (x=0, 플라즈마 열유속 경계조건)
        double qSurf = plasmaHeatFluxWattsM2;
        double tSurfOld = _temperatures[0];
        _temperatures[0] = tSurfOld + 2.0 * r * (_temperatures[1] - tSurfOld + (qSurf * Dx / ThermalCondW));

        // 내부 노드 전도 계산
        for (int i = 1; i < GridNodes - 1; i++)
        {
            double tPrev = _temperatures[i - 1];
            double tCurr = _temperatures[i];
            double tNext = _temperatures[i + 1];
            _temperatures[i] = tCurr + r * (tPrev - 2.0 * tCurr + tNext);
        }

        // 냉각 채널 노드 (x=L, 대류 경계조건)
        int last = GridNodes - 1;
        double tLastOld = _temperatures[last];
        _temperatures[last] = tLastOld + 2.0 * r * (_temperatures[last - 1] - tLastOld - (hCool * Dx / ThermalCondW) * (tLastOld - _coolantTemp));

        // 3. 열응력 계산
        double tSurf = _temperatures[0];
        double tAvg = 0.0;
        for (int i = 0; i < GridNodes; i++)
        {
            tAvg += _temperatures[i];
        }
        tAvg /= GridNodes;

        double deltaT = tSurf - tAvg;
        double thermalStress = (ElasticModulus * ThermalExpansion / (1.0 - PoissonRatio)) * deltaT;
        double thermalStressMpa = thermalStress * 1e-6;

        // 4. 냉각재 유량 능동 피드백 제어 (표면 1,000°C 초과 시 유량 급증)
        double targetFlow = 5.0;
        if (tSurf > 900.0)
        {
            double excess = tSurf - 900.0;
            targetFlow += excess * 0.05; // 비례 냉각 증강
        }
        if (targetFlow > 25.0) targetFlow = 25.0; // 펌프 최대 용량 제한
        _currentCoolantFlow = _currentCoolantFlow * 0.9 + targetFlow * 0.1; // 1차 필터 응답

        bool risk = tSurf >= RecrystallizationLimit;

        return new DivertorStatus(tSurf, _temperatures[last], thermalStressMpa, _currentCoolantFlow, risk);
    }
}