BWRSimlator/Reactor.cs

namespace BWR_siumulator
{
    public class ThermalPort
    {
        public double Measurement { get; set; }
        public double ValueToAffect { get; set; }

        public ThermalPort(double measurement, double valueToAffect)
        {
            Measurement = measurement;
            ValueToAffect = valueToAffect;
        }
    }

    public class Reactor
    {
        // --- Core Physics and Material Constants ---
        public const double CONTROL_ROD_SPEED = 1.0;
        public const double DECAY_HEAT_FRACTION = 0.065;
        public const double DECAY_HEAT_LAMBDA = 0.0077;
        public const double BETA = 0.0067;
        public const double LAMBDA_DECAY = 0.078;
        public const double L_PROMPT_NEUTRON = 0.00002;
        public const double FUEL_HEAT_CAPACITY = 2.8;
        public const double HEAT_TRANSFER_COEFF = 0.5;
        public const double ALPHA_SQRT_T_FUEL = -0.008;
        public const double CONTROL_ROD_WORTH = -0.05;
        public const double INTRINSIC_SOURCE = 0.25;
        public const double COOLANT_TEMP = 300; // K
        public const double STABLE_TEMP = 547; // K [The temperature at which the doppler effect has no effect on reactivity]

        // --- State Variables ---
        public double Power { get; private set; } // MW

        // --- Inspect Variables ---
        public double ReactivityDoppler { get; private set; }
        public double ReactivityRod { get; private set; }
        
        public double FuelTemp { get; private set; } // K
        public double HeatGenerated { get; private set; }
        
        public double Precursors { get; private set; }
        public double DecayHeatPrecursors { get; private set; }
        
        public double ReactivityTotal { get; private set; }
        public double ControlRodCurrentPosition { get; private set; }
        public double ControlRodTargetPosition { get; set; } // Can be set externally
        public double ReactorPeriod { get; private set; }

        public ThermalPort ThermalPort { get; private set; }

        public Reactor(double initialPower = 1.0, double fuelTemp = 300)
        {
            Power = initialPower;
            FuelTemp = fuelTemp;

            // Initialize precursors and decay heat to be in equilibrium
            Precursors = Power * BETA / (LAMBDA_DECAY * L_PROMPT_NEUTRON);
            DecayHeatPrecursors = Power * DECAY_HEAT_FRACTION / DECAY_HEAT_LAMBDA;
            
            ReactivityTotal = 0.0;
            ControlRodCurrentPosition = 100;
            ControlRodTargetPosition = 100;
            ReactorPeriod = double.PositiveInfinity;

            ThermalPort = new ThermalPort(FuelTemp, HeatGenerated);
        }
        private double GetReactivityFromRodPosition()
        {
            double positionRad = (ControlRodCurrentPosition / 100.0) * Math.PI;
            double effectiveness = (1 - Math.Cos(positionRad)) / 2.0;
            return CONTROL_ROD_WORTH * effectiveness;
        }

        /// <summary>
        /// Calculates reactivity from the fuel temperature.
        /// </summary>
        /// <returns>The reactivity due to fuel temperature (Doppler effect).</returns>
        private double CalculateDopplerReactivity()
        {
            return ALPHA_SQRT_T_FUEL * (Math.Sqrt(FuelTemp) - Math.Sqrt(STABLE_TEMP));
        }

        /// <summary>
        /// Advances the simulation by one time step, dt.
        /// </summary>
        /// <param name="dt">The time step in seconds.</param>
        public void Step(double dt)
        {
            // --- UPDATE CONTROL ROD POSITION ---
            if (ControlRodCurrentPosition != ControlRodTargetPosition)
            {
                double difference = ControlRodTargetPosition - ControlRodCurrentPosition;
                double maxMove = CONTROL_ROD_SPEED * dt;
                if (Math.Abs(difference) < maxMove)
                {
                    ControlRodCurrentPosition = ControlRodTargetPosition;
                }
                else if (difference > 0)
                {
                    ControlRodCurrentPosition += maxMove;
                }
                else
                {
                    ControlRodCurrentPosition -= maxMove;
                }
            }

            // --- PHYSICS CALCULATION ---
            // 1. Calculate Total Reactivity
            ReactivityDoppler = CalculateDopplerReactivity();
            ReactivityRod = GetReactivityFromRodPosition();
            ReactivityTotal = ReactivityDoppler + ReactivityRod;

            // 2. Solve Point Kinetics
            double powerOld = Power;
            double dtLambda = dt * LAMBDA_DECAY;
            double dtOverL = dt / L_PROMPT_NEUTRON;
            double numerator = Power + Precursors * dtLambda / (1 + dtLambda) + dt * INTRINSIC_SOURCE;
            double denominator = 1 - dtOverL * (ReactivityTotal - BETA) -
                                 (dtOverL * BETA * dtLambda) / (1 + dtLambda);
            Power = numerator / denominator;
            
            Precursors = (Precursors + dt * Power * BETA / L_PROMPT_NEUTRON) / (1 + dtLambda);
            
            // 3. Calculate Reactor Period
            double powerChange = Power - powerOld;
            if (Math.Abs(powerChange) > 1e-9 && Power > 1e-7)
            {
                ReactorPeriod = (Power * dt) / powerChange;
                if (ReactorPeriod < -1500 || ReactorPeriod > 1500)
                {
                    ReactorPeriod = double.PositiveInfinity;
                }
            }
            else
            {
                ReactorPeriod = double.PositiveInfinity;
            }

            // 4. Solve for Decay Heat
            DecayHeatPrecursors = (DecayHeatPrecursors + dt * Power * DECAY_HEAT_FRACTION) / (1 + dt * DECAY_HEAT_LAMBDA);

            // 5. Solve Thermal-Hydraulics
            double promptHeat = Power * (1 - DECAY_HEAT_FRACTION);
            HeatGenerated = promptHeat + DecayHeatPrecursors;
            //double heatRemoved = HEAT_TRANSFER_COEFF * (FuelTemp - COOLANT_TEMP);
            //double dFuelTemp = (HeatGenerated - heatRemoved) / FUEL_HEAT_CAPACITY;
            double dFuelTemp = HeatGenerated / FUEL_HEAT_CAPACITY;
            FuelTemp += dFuelTemp * dt;
            
            // --- CLAMPING ---
            if (Power < 0) Power = 0;
            if (Precursors < 0) Precursors = 0;
            if (DecayHeatPrecursors < 0) DecayHeatPrecursors = 0;

            // Update ThermalPort
            ThermalPort.Measurement = FuelTemp;
            ThermalPort.ValueToAffect = HeatGenerated;
        }
    }
}