도산안창호급(KSS-III) 잠수함발사탄도미사일(SLBM) 수중 캐니스터 콜드런치 캡슐 이탈 기포막(Supercavitation Bubble Screen) 유동 & 해수면 수직 돌파(Broaching) 6자유도 1,000Hz 0-Alloc C# 엔진

using System;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;

namespace ZeliDesk.Defense.SubmarineWeapons
{
    public enum SlbmPhase
    {
        InTubeGasExpansion,     // 캐니스터 내 가스 팽창 사출
        UnderwaterSupercavitating, // 수중 초공동 기포막 상승
        SurfaceBroaching,       // 해수면 돌파 슬래밍 구간
        AirborneAttitudeControl,// 공중 수직 자세 안정화
        MainMotorIgnited        // 주 추진기관 점화 상승
    }

    [StructLayout(LayoutKind.Sequential, Pack = 16)]
    public struct SlbmState
    {
        public Vector256<float> PosVel;       // [x, y, z, 0, vx, vy, vz, 0] (z 음수: 수심 m, vz: 수직속도 m/s)
        public Vector256<float> AnglesRate;   // [roll, pitch, yaw, 0, p, q, r, 0]
        public float CavitationNumber;        // 캐비테이션 수 sigma
        public float AmbientPressureKPa;      // 주변 압력 (kPa)
        public float HydrodynamicDragN;       // 유체 저항력 (N)
        public SlbmPhase Phase;               // 현재 비행 단계
    }

    public sealed class Kss3SlbmLaunchEngine
    {
        private const float MissileMass = 18500.0f;     // 미사일 총질량 (kg)
        private const float Diameter = 1.4f;            // 동체 직경 (m)
        private const float Length = 9.8f;              // 전장 (m)
        private const float RhoWater = 1025.0f;         // 해수 밀도 (kg/m^3)
        private const float RhoAir = 1.225f;            // 대기 밀도 (kg/m^3)
        private const float Gravity = 9.80665f;
        private const float AreaBase = MathF.PI * (Diameter * 0.5f) * (Diameter * 0.5f);

        [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
        public void StepSlbmLaunch1000Hz(ref SlbmState state, float dt)
        {
            float z = state.PosVel.GetElement(2); // 수심 (m, 수면: 0, 수중: < 0)
            float vz = state.PosVel.GetElement(6); // 상승 속도 (m/s)
            float pitch = state.AnglesRate.GetElement(1);
            float q = state.AnglesRate.GetElement(5);

            // 1. 수심에 따른 주변 정수압 계산 (P = P_atm + rho * g * h)
            float depthM = MathF.Max(-z, 0.0f);
            float hydrostaticKPa = 101.325f + (RhoWater * Gravity * depthM) * 0.001f;
            state.AmbientPressureKPa = hydrostaticKPa;

            // 2. 비행 단계 판정 (FSM)
            if (depthM > 35.0f)
            {
                state.Phase = SlbmPhase.InTubeGasExpansion;
            }
            else if (depthM > 2.0f)
            {
                state.Phase = SlbmPhase.UnderwaterSupercavitating;
            }
            else if (depthM >= 0.0f || z < 3.0f)
            {
                state.Phase = SlbmPhase.SurfaceBroaching;
            }
            else if (z >= 3.0f && z < 18.0f)
            {
                state.Phase = SlbmPhase.AirborneAttitudeControl;
            }
            else
            {
                state.Phase = SlbmPhase.MainMotorIgnited;
            }

            // 3. 단계별 유체력 및 추력 계산
            float netForceZ = 0f;
            float pitchTorque = 0f;

            switch (state.Phase)
            {
                case SlbmPhase.InTubeGasExpansion:
                    // 사출 가스 압축 팽창력 (초기 2.5 MPa 가스압)
                    float gasPressureKPa = 2400.0f * MathF.Max(1.0f - (45.0f - depthM) / 10.0f, 0.3f);
                    float ejectForce = (gasPressureKPa - hydrostaticKPa) * 1000.0f * AreaBase;
                    netForceZ = ejectForce - MissileMass * Gravity;
                    break;

                case SlbmPhase.UnderwaterSupercavitating:
                    // 초공동 기포막 형성 (캐비테이션 수 sigma ~ 0.08)
                    float dynamicPressure = 0.5f * RhoWater * vz * vz;
                    state.CavitationNumber = (hydrostaticKPa - 35.0f) * 1000.0f / MathF.Max(dynamicPressure, 1.0f);

                    // 기포막으로 인한 마찰 항력 85% 감쇄
                    float cdSupercav = 0.065f;
                    float dragWater = cdSupercav * dynamicPressure * AreaBase;
                    state.HydrodynamicDragN = dragWater;

                    // 아르키메데스 부력
                    float displacedVol = MathF.PI * (Diameter * 0.5f) * (Diameter * 0.5f) * Length * 0.75f;
                    float buoyancy = RhoWater * Gravity * displacedVol;

                    netForceZ = buoyancy - MissileMass * Gravity - dragWater;
                    break;

                case SlbmPhase.SurfaceBroaching:
                    // 해수면 통과 시 유체 충격 슬래밍(Slamming) 및 비대칭 전복 모멘트
                    float wettedFraction = Math.Clamp(-z / Length, 0.0f, 1.0f);
                    float mixedDensity = wettedFraction * RhoWater + (1.0f - wettedFraction) * RhoAir;
                    float broachDrag = 0.25f * (0.5f * mixedDensity * vz * vz) * AreaBase;
                    state.HydrodynamicDragN = broachDrag;

                    // 파랑에 의한 외란 피칭 토크 (파고 2.5m 모델)
                    pitchTorque = 85000.0f * MathF.Sin(depthM * 2.0f) - 15000.0f * q;
                    netForceZ = -MissileMass * Gravity - broachDrag;
                    break;

                case SlbmPhase.AirborneAttitudeControl:
                    // 공중 자세 안정화 (보조 펄스 측추력기 작동)
                    state.HydrodynamicDragN = 0.5f * RhoAir * vz * vz * AreaBase * 0.15f;
                    // 피치 각도 수직 복원 PD 제어
                    pitchTorque = -120000.0f * pitch - 25000.0f * q;
                    netForceZ = -MissileMass * Gravity;
                    break;

                case SlbmPhase.MainMotorIgnited:
                    // 1단 주 로켓 모터 점화 (추력 420 kN)
                    float rocketThrust = 420000.0f;
                    netForceZ = rocketThrust - MissileMass * Gravity;
                    pitchTorque = -80000.0f * pitch - 20000.0f * q;
                    break;
            }

            // 4. 상태 적분 (RK2)
            float accelZ = netForceZ / MissileMass;
            float qDot = pitchTorque / 185000.0f; // 피치 관성모멘트 (kg*m^2)

            float newVz = vz + accelZ * dt;
            float newZ = z + 0.5f * (vz + newVz) * dt;

            float newQ = q + qDot * dt;
            float newPitch = pitch + 0.5f * (q + newQ) * dt;

            state.PosVel = Vector256.Create(
                state.PosVel.GetElement(0),
                state.PosVel.GetElement(1),
                newZ, 0f,
                state.PosVel.GetElement(4),
                state.PosVel.GetElement(5),
                newVz, 0f);

            state.AnglesRate = Vector256.Create(
                state.AnglesRate.GetElement(0),
                newPitch,
                state.AnglesRate.GetElement(2),
                0f,
                state.AnglesRate.GetElement(4),
                newQ,
                state.AnglesRate.GetElement(6),
                0f);
        }
    }
}