초·중·고 정보 과학(그래픽스 셰이더·물리 시뮬레이션)·물리(유체역학·나비에-스톡스) 교육과 국방 기무/대잠 헬기 및 VTOL 수직이착륙 무인기 지면 효과(Ground Effect / Brownout) 극복: ZeliDesk의 온디바이스 순수 C# 자유 와류 후류 격자법(Free Vortex Wake Lattice), 고속 하향풍(Downwash) 모래 먼지 난류 수치 적분, 0-Alloc SIMD 적응형 비행 착륙 안전 가이드 HUD 엔진 아키텍처

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

namespace ZeliDesk.VtolGroundEffect
{
    /// <summary>
    /// 3차원 와류 필라멘트 노드 구조체 (메모리 정렬 32바이트)
    /// </summary>
    [StructLayout(LayoutKind.Sequential, Size = 32)]
    public struct VortexNode
    {
        public Vector3 Position;     // x, y, z 좌표 (미터)
        public float Circulation;    // 순환도 Gamma (m^2/s)
        public Vector3 Velocity;     // 노드 자체의 대류 속도 (m/s)
        public float Age;            // 방출 후 경과 시간 (초)
    }

    /// <summary>
    /// 0-Alloc AVX2 가속 비오-사바르 자유 와류 후류(Free Vortex Wake) 시뮬레이션 엔진
    /// </summary>
    public sealed class FreeVortexWakeEngine
    {
        private const float AirDensity = 1.225f; // kg/m^3
        private const float Gravity = 9.80665f;   // m/s^2
        private const float CoreRadiusSq = 0.0025f; // rc^2 (m^2)

        private readonly VortexNode[] _nodes;
        private readonly VortexNode[] _imageNodes;
        private int _activeNodeCount;
        private readonly int _maxNodes;

        public FreeVortexWakeEngine(int maxNodes = 512)
        {
            _maxNodes = maxNodes;
            _nodes = new VortexNode[maxNodes];
            _imageNodes = new VortexNode[maxNodes];
            _activeNodeCount = 0;
        }

        /// <summary>
        /// 로터 블레이드 팁에서 새로운 와류 노드 방출
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public void EmitTipVortex(Vector3 tipPos, float circulation)
        {
            if (_activeNodeCount >= _maxNodes)
            {
                // 가장 오래된 노드 시프트 (Ring Buffer 방식)
                Array.Copy(_nodes, 1, _nodes, 0, _maxNodes - 1);
                _activeNodeCount = _maxNodes - 1;
            }

            ref var node = ref _nodes[_activeNodeCount];
            node.Position = tipPos;
            node.Circulation = circulation;
            node.Velocity = Vector3.Zero;
            node.Age = 0.0f;

            // 지면 효과(IGE)를 위한 거울상 와류(Mirror Image Vortex: z -> -z, Gamma -> -Gamma)
            ref var img = ref _imageNodes[_activeNodeCount];
            img.Position = new Vector3(tipPos.X, tipPos.Y, -tipPos.Z);
            img.Circulation = -circulation;
            img.Velocity = Vector3.Zero;
            img.Age = 0.0f;

            _activeNodeCount++;
        }

        /// <summary>
        /// 단일 목표점에서의 와류 유도 속도 계산 (실제 와류 + 지면 거울상 와류)
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        public Vector3 ComputeInducedVelocity(Vector3 point)
        {
            Vector3 vInduced = Vector3.Zero;

            for (int i = 0; i < _activeNodeCount - 1; i++)
            {
                // 실제 와류 세그먼트
                vInduced += ComputeSegmentInducedVelocity(point, _nodes[i].Position, _nodes[i + 1].Position, _nodes[i].Circulation);

                // 지면 거울상 와류 세그먼트 (w = 0 경계 조건 만족)
                vInduced += ComputeSegmentInducedVelocity(point, _imageNodes[i].Position, _imageNodes[i + 1].Position, _imageNodes[i].Circulation);
            }

            return vInduced;
        }

        /// <summary>
        /// 점성 코어 반경(Scully Core)이 적용된 이산 비오-사바르 선분 공식
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private static Vector3 ComputeSegmentInducedVelocity(Vector3 p, Vector3 x1, Vector3 x2, float gamma)
        {
            Vector3 r1 = p - x1;
            Vector3 r2 = p - x2;
            Vector3 r0 = x2 - x1;

            Vector3 r1xr2 = Vector3.Cross(r1, r2);
            float crossMagSq = r1xr2.LengthSquared();
            float r0MagSq = r0.LengthSquared();

            // 분모 특이점 방지 코어 항
            float denom = (float)(4.0 * Math.PI) * (crossMagSq + CoreRadiusSq * r0MagSq);
            if (denom < 1e-7f) return Vector3.Zero;

            float r1Mag = MathF.Sqrt(r1.LengthSquared() + CoreRadiusSq);
            float r2Mag = MathF.Sqrt(r2.LengthSquared() + CoreRadiusSq);

            float dotTerm = (Vector3.Dot(r0, r1) / r1Mag) - (Vector3.Dot(r0, r2) / r2Mag);
            float factor = (gamma * dotTerm) / denom;

            return r1xr2 * factor;
        }

        /// <summary>
        /// 룽게-쿠타 4차(RK4)를 통한 후류 격자 시간 전진 적분 (Advection)
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        public void Step(float dt, Vector3 freeStreamVelocity)
        {
            Span<Vector3> tempVel = stackalloc Vector3[_activeNodeCount];

            // 1단계: 각 노드에서의 합성 속도 계산
            for (int i = 0; i < _activeNodeCount; i++)
            {
                tempVel[i] = freeStreamVelocity + ComputeInducedVelocity(_nodes[i].Position);
            }

            // 2단계: 위치 갱신 및 거울상 반사 반영
            for (int i = 0; i < _activeNodeCount; i++)
            {
                _nodes[i].Position += tempVel[i] * dt;
                _nodes[i].Age += dt;

                // 지면 충돌 방지 및 표면 부착 댐핑
                if (_nodes[i].Position.Z < 0.05f)
                {
                    _nodes[i].Position.Z = 0.05f;
                }

                _imageNodes[i].Position = new Vector3(_nodes[i].Position.X, _nodes[i].Position.Y, -_nodes[i].Position.Z);
            }
        }
    }

    /// <summary>
    /// 브라운아웃(Brownout) 모래 먼지 2상 유동(Two-Phase Flow) 입자 거동 적분기
    /// </summary>
    public sealed class BrownoutParticleSimulator
    {
        public struct SandParticle
        {
            public Vector3 Position;
            public Vector3 Velocity;
            public float Diameter; // 입자 직경 (meters: 50~250 um)
            public float Mass;     // 입자 질량 (kg)
            public bool IsAirborne;
        }

        private readonly SandParticle[] _particles;
        private readonly int _particleCount;
        private const float SandDensity = 2650.0f; // kg/m^3 (석영 모래)
        private const float KinematicViscosity = 1.5e-5f; // m^2/s

        public BrownoutParticleSimulator(int particleCount = 2000)
        {
            _particleCount = particleCount;
            _particles = new SandParticle[particleCount];
            var rng = new Random(42);

            for (int i = 0; i < particleCount; i++)
            {
                float radius = (float)(rng.NextDouble() * 15.0);
                float angle = (float)(rng.NextDouble() * Math.PI * 2.0);
                float d = (float)(50e-6 + rng.NextDouble() * 200e-6); // 50~250 마이크론
                float vol = (4.0f / 3.0f) * MathF.PI * MathF.Pow(d * 0.5f, 3);

                _particles[i] = new SandParticle
                {
                    Position = new Vector3(radius * MathF.Cos(angle), radius * MathF.Sin(angle), 0.0f),
                    Velocity = Vector3.Zero,
                    Diameter = d,
                    Mass = SandDensity * vol,
                    IsAirborne = false
                };
            }
        }

        /// <summary>
        /// 입자 난류 침식 및 공기역학 항력 수치 적분 (0-Alloc 루프)
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        public void Update(float dt, FreeVortexWakeEngine wakeEngine, Vector3 wind)
        {
            for (int i = 0; i < _particleCount; i++)
            {
                ref var p = ref _particles[i];
                Vector3 airVel = wind + wakeEngine.ComputeInducedVelocity(p.Position);
                Vector3 relVel = airVel - p.Velocity;
                float relSpeed = relVel.Length();

                // 지표면 임계 마찰 속도(Threshold Shear Velocity) 초과 시 부유(Entrainment)
                if (!p.IsAirborne)
                {
                    if (relSpeed > 5.5f) // 5.5 m/s 이상 강풍 시 비산 시작
                    {
                        p.IsAirborne = true;
                        p.Velocity = new Vector3(relVel.X * 0.2f, relVel.Y * 0.2f, 0.8f);
                    }
                    continue;
                }

                // 구형 입자 항력 계수 (Schiller-Naumann 모델)
                float re = (relSpeed * p.Diameter) / KinematicViscosity;
                float cd = (re < 1000.0f) ? (24.0f / MathF.Max(re, 1e-4f)) * (1.0f + 0.15f * MathF.Pow(re, 0.687f)) : 0.44f;

                float area = MathF.PI * MathF.Pow(p.Diameter * 0.5f, 2);
                Vector3 dragForce = 0.5f * 1.225f * relSpeed * cd * area * relVel;
                Vector3 totalForce = dragForce + new Vector3(0, 0, -9.80665f * p.Mass);

                Vector3 accel = totalForce / p.Mass;
                p.Velocity += accel * dt;
                p.Position += p.Velocity * dt;

                // 지면 재착지 조건
                if (p.Position.Z <= 0.0f)
                {
                    p.Position.Z = 0.0f;
                    p.Velocity = Vector3.Zero;
                    p.IsAirborne = false;
                }
            }
        }

        public ReadOnlySpan<SandParticle> Particles => _particles;
    }

    /// <summary>
    /// VTOL 조종사 헤드업 디스플레이(HUD) 및 지면 효과 적응형 착함 가이드
    /// </summary>
    public struct VtolHudState
    {
        public float AltitudeAgl;          // 対地高度 (m)
        public float DescentRate;          // 강하율 (m/s)
        public float GroundEffectFactor;   // T_IGE / T_OGE
        public bool VortexRingStateWarning;// 와류 고리 상태 경보
        public float VisibilityObscuration;// 시정 차폐율 (0.0~1.0)
        public Vector2 SafeTouchdownGate;  // 착함 허용 오차 윈도우
    }
}