FluidSim/Scenarios/Simulation.cs

using FluidSim.Audio;
using FluidSim.Components;
using FluidSim.Core;
using FluidSim.Utilities;

namespace FluidSim.Scenarios
{
    public static class Simulation
    {
        private static Solver solver;
        private static Volume0D cavity, ambient;
        private static Pipe1D neck;
        private static Connection connNeckCavity, connNeckAmbient;
        private static int stepCount;
        private static double time;
        private static double dt;

        public static void Initialize(int sampleRate)
        {
            dt = 1.0 / sampleRate;

            // --- Cavity (the “bottle”) ---
            double V_cav = 1.0 * Units.L;                 // 1 litre
            cavity = new Volume0D(V_cav, 10 * Units.atm, Units.Celsius(20), sampleRate);

            // --- Ambient (a huge “room”) ---
            double V_amb = 1000.0;                        // 1000 m³ ≈ constant pressure
            ambient = new Volume0D(V_amb, 1.0 * Units.atm, Units.Celsius(20), sampleRate);

            // --- Neck (short pipe) ---
            double length_neck = 80.0 * Units.mm;
            double diam_neck = 10 * Units.mm;
            double area_neck = Units.AreaFromRadius(diam_neck, Units.mm);                             // few cells enough for a short neck
            neck = new Pipe1D(length_neck, area_neck, sampleRate);
            neck.SetUniformState(ambient.Density, 0.0, ambient.Pressure);
            neck.FrictionFactor = 0.02;                    // slight damping

            // --- Connections ---
            // Cavity-to-neck (full area)
            connNeckCavity = new Connection(cavity.Port, neck.PortA)
            {
                Area = area_neck,
                DischargeCoefficient = 1.0,
                Gamma = 1.4
            };
            // Neck-to-ambient (SoundConnection to capture the radiated tone)
            connNeckAmbient = new SoundConnection(neck.PortB, ambient.Port)
            {
                Area = area_neck,
                DischargeCoefficient = 1.0,
                Gamma = 1.4
            };

            // --- Solver ---
            solver = new Solver();
            solver.AddVolume(cavity);
            solver.AddVolume(ambient);
            solver.AddPipe(neck);
            solver.AddConnection(connNeckCavity);
            solver.AddConnection(connNeckAmbient);

            // --- Sound tuning ---
            SoundProcessor.MaxDeltaP = 0.1f * (float)Units.atm; // small Δp expected
            SoundProcessor.MaxArea = (float)area_neck;
            SoundProcessor.MaxVelocity = 343f;
            SoundProcessor.ReferenceDensity = 1.2f;
            SoundProcessor.ReferenceSpeedOfSound = 343f;
            SoundProcessor.Gain = 10.0f;  // amplify because Δp is small
        }

        public static float Process()
        {
            solver.Step();
            time += dt;
            stepCount++;
            Log();
            return solver.LastSample;
        }

        public static void Log()
        {
            if (stepCount <= 200 || stepCount % (int)(0.5 / dt) == 0)
            {
                Console.WriteLine(
                    $"t = {time * 1e3:F3} ms  " +
                    $"Sample = {solver.LastSample:F6},  " +
                    $"P_cav = {cavity.Pressure / 1e5:F6} bar,  " +
                    $"flow_cav = {cavity.Port.MassFlowRate / 1e3:F4} g/s,  " +
                    $"flow_neck = {neck.PortB.MassFlowRate * 1e3:F4} g/s");
            }
        }
    }
}