twine-models 0.4.0

//! Isentropic compressor model.
//!
//! Computes an outlet state at a target pressure using an isentropic efficiency `eta`.
//!
//! Given an inlet state and a target outlet pressure, the model:
//! 1. computes an ideal (isentropic) outlet at constant entropy: `(p_out, s_in)`,
//! 2. evaluates the ideal enthalpy rise `Δh_s = h_out,s − h_in`,
//! 3. maps to an actual enthalpy rise using `η = Δh_s / Δh` ⇒ `Δh = Δh_s / η`,
//! 4. requests an outlet state from `(p_out, h_out)` and reports the required work.

use uom::si::f64::Pressure;

use crate::support::{
    thermo::{
        State,
        capability::{HasEnthalpy, HasEntropy, HasPressure, StateFrom, ThermoModel},
    },
    turbomachinery::{
        CompressionWork, InletProperties, IsentropicEfficiency,
        compressor::{CompressionError, CompressionResult},
    },
    units::{SpecificEnthalpy, SpecificEntropy},
};

/// Computes the compressor outlet state and required work using an isentropic efficiency model.
///
/// `eta` is an isentropic efficiency in `(0, 1]`.
/// Values near zero represent extremely inefficient compression and result in
/// very large work.
///
/// The returned [`CompressionResult::work`] is the target enthalpy rise
/// `h_out − h_in`, reported as a [`CompressionWork`]. If the thermodynamic
/// model returns the outlet state via inverse-solving or interpolation,
/// `thermo.enthalpy(&outlet)` may differ slightly from the target.
///
/// # Errors
///
/// Returns [`CompressionError`] if the thermodynamic model fails, `p_out < p_in`,
/// or the resulting work is non-physical.
pub fn isentropic<Fluid, Model>(
    inlet: &State<Fluid>,
    p_out: Pressure,
    eta: IsentropicEfficiency,
    thermo: &Model,
) -> Result<CompressionResult<Fluid>, CompressionError<Fluid>>
where
    Fluid: Clone,
    Model: ThermoModel<Fluid = Fluid>
        + HasPressure
        + HasEnthalpy
        + HasEntropy
        + StateFrom<(Fluid, Pressure, SpecificEnthalpy)>
        + StateFrom<(Fluid, Pressure, SpecificEntropy)>,
{
    let p_in = thermo
        .pressure(inlet)
        .map_err(CompressionError::inlet_pressure_failed)?;

    let h_in = thermo
        .enthalpy(inlet)
        .map_err(CompressionError::inlet_enthalpy_failed)?;

    let s_in = thermo
        .entropy(inlet)
        .map_err(CompressionError::inlet_entropy_failed)?;

    let inlet_props = InletProperties {
        thermo,
        fluid: inlet.fluid.clone(),
        p_in,
        h_in,
        s_in,
    };
    isentropic_core(inlet_props, p_out, eta)
}

/// Core isentropic compression model.
///
/// # Errors
///
/// Returns [`CompressionError`] if the thermodynamic model fails, `p_out < p_in`,
/// or the resulting work is non-physical.
#[allow(clippy::needless_pass_by_value)]
fn isentropic_core<Fluid, Model>(
    inlet_props: InletProperties<'_, Fluid, Model>,
    p_out: Pressure,
    eta: IsentropicEfficiency,
) -> Result<CompressionResult<Fluid>, CompressionError<Fluid>>
where
    Model: ThermoModel<Fluid = Fluid>
        + HasEnthalpy
        + StateFrom<(Fluid, Pressure, SpecificEnthalpy)>
        + StateFrom<(Fluid, Pressure, SpecificEntropy)>,
{
    let InletProperties {
        thermo,
        fluid,
        p_in,
        h_in,
        s_in,
    } = inlet_props;

    if p_out < p_in {
        return Err(CompressionError::OutletPressureLessThanInlet { p_in, p_out });
    }

    let isentropic_outlet = thermo.state_from((fluid, p_out, s_in)).map_err(|source| {
        CompressionError::ideal_outlet_state_from_pressure_entropy_failed(p_out, s_in, source)
    })?;

    let h_out_s = thermo
        .enthalpy(&isentropic_outlet)
        .map_err(CompressionError::ideal_outlet_enthalpy_failed)?;

    let State { fluid, .. } = isentropic_outlet;

    let dh_s = h_out_s - h_in;
    let dh_actual = dh_s / eta.ratio();
    let h_out_target = h_in + dh_actual;

    let outlet = thermo
        .state_from((fluid, p_out, h_out_target))
        .map_err(|source| {
            CompressionError::outlet_state_from_pressure_enthalpy_failed(
                p_out,
                h_out_target,
                source,
            )
        })?;

    let raw_work = h_out_target - h_in;

    match CompressionWork::new(raw_work) {
        Ok(work) => Ok(CompressionResult { outlet, work }),
        Err(_) => Err(CompressionError::NonPhysicalWork { outlet, raw_work }),
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    use approx::assert_relative_eq;
    use uom::si::{
        f64::{MassDensity, Pressure, ThermodynamicTemperature},
        mass_density::kilogram_per_cubic_meter,
        pressure::{kilopascal, pascal},
        specific_heat_capacity::joule_per_kilogram_kelvin,
        thermodynamic_temperature::kelvin,
    };

    use crate::support::{
        thermo::{PropertyError, capability::HasPressure},
        turbomachinery::{
            IsentropicEfficiency,
            test_utils::{FakeMode, FakeThermo, MockGas, enth_si, mock_gas_model},
        },
        units::SpecificEntropy,
    };

    #[test]
    fn normal_compression_matches_expected_work() {
        let thermo = mock_gas_model();

        let p_in = Pressure::new::<kilopascal>(100.0);
        let t_in = ThermodynamicTemperature::new::<kelvin>(300.0);
        let inlet: State<MockGas> = thermo.state_from((MockGas, t_in, p_in)).unwrap();

        // Pick `p2/p1 = 2^7` so that `T2s = T1*(p2/p1)^((k-1)/k) = T1*(2^7)^(2/7) = 4*T1`.
        let p_out = Pressure::new::<kilopascal>(12_800.0);
        let eta = IsentropicEfficiency::new(0.9).unwrap();

        let result = isentropic(&inlet, p_out, eta, &thermo).unwrap();

        // Expected output for an ideal gas:
        // ```
        // k = cp/(cp-R) = 1.4
        // T2s = 4*T1 = 1200 K
        // w = cp*(T2s - T1)/eta = 1_000_000 J/kg
        // ```
        let expected_work_j_per_kg = 1_000_000.0;
        assert_relative_eq!(result.work.quantity().value, expected_work_j_per_kg);

        // The outlet state is computed from (p_out, h_out_target),
        // so for an ideal gas: `T2 = T1 + w/cp = 300 K + 1_000_000/1000 = 1300 K`.
        assert_relative_eq!(result.outlet.temperature.get::<kelvin>(), 1300.0);
    }

    #[test]
    fn outlet_pressure_less_than_inlet_is_an_error() {
        let thermo = mock_gas_model();

        let inlet = State {
            temperature: ThermodynamicTemperature::new::<kelvin>(300.0),
            density: MassDensity::new::<kilogram_per_cubic_meter>(1.0),
            fluid: MockGas,
        };
        let p_in = thermo.pressure(&inlet).unwrap();
        let p_out = p_in - Pressure::new::<pascal>(1.0);

        let err = isentropic(
            &inlet,
            p_out,
            IsentropicEfficiency::new(1.0).unwrap(),
            &thermo,
        )
        .unwrap_err();

        match err {
            CompressionError::OutletPressureLessThanInlet {
                p_in: reported_p_in,
                p_out: reported_p_out,
            } => {
                assert_eq!(reported_p_in, p_in);
                assert_eq!(reported_p_out, p_out);
            }
            other => panic!("unexpected error: {other:?}"),
        }
    }

    // Shared inputs for `isentropic_core` error-path tests.
    fn core_inputs() -> (
        Pressure,
        SpecificEnthalpy,
        SpecificEntropy,
        Pressure,
        IsentropicEfficiency,
    ) {
        let p_in = Pressure::new::<kilopascal>(100.0);
        let h_in = enth_si(0.0);
        let s_in = SpecificEntropy::new::<joule_per_kilogram_kelvin>(0.0);
        let p_out = Pressure::new::<kilopascal>(200.0);
        let eta = IsentropicEfficiency::new(1.0).unwrap();
        (p_in, h_in, s_in, p_out, eta)
    }

    #[test]
    fn core_state_from_pressure_entropy_failure_is_wrapped() {
        let thermo = FakeThermo {
            mode: FakeMode::FailStateFromPressureEntropy,
        };
        let (p_in, h_in, s_in, p_out, eta) = core_inputs();

        let inlet_props = InletProperties {
            thermo: &thermo,
            fluid: MockGas,
            p_in,
            h_in,
            s_in,
        };
        let err = isentropic_core(inlet_props, p_out, eta).unwrap_err();

        match err {
            CompressionError::ThermodynamicModelFailed { context, source: _ } => {
                assert!(context.contains("ideal_outlet_state_from("));
                assert!(context.contains("p_out="));
                assert!(context.contains("s_in="));
            }
            other => panic!("unexpected error: {other:?}"),
        }
    }

    #[test]
    fn core_state_from_pressure_enthalpy_failure_is_wrapped() {
        let thermo = FakeThermo {
            mode: FakeMode::FailStateFromPressureEnthalpy,
        };
        let (p_in, h_in, s_in, p_out, eta) = core_inputs();

        let inlet_props = InletProperties {
            thermo: &thermo,
            fluid: MockGas,
            p_in,
            h_in,
            s_in,
        };
        let err = isentropic_core(inlet_props, p_out, eta).unwrap_err();

        match err {
            CompressionError::ThermodynamicModelFailed { context, source: _ } => {
                assert!(context.contains("outlet_state_from("));
                assert!(context.contains("p_out="));
                assert!(context.contains("h_out_target="));
            }
            other => panic!("unexpected error: {other:?}"),
        }
    }

    #[test]
    fn core_property_failure_is_wrapped() {
        let thermo = FakeThermo {
            mode: FakeMode::FailEnthalpy,
        };
        let (p_in, h_in, s_in, p_out, eta) = core_inputs();

        let inlet_props = InletProperties {
            thermo: &thermo,
            fluid: MockGas,
            p_in,
            h_in,
            s_in,
        };
        let err = isentropic_core(inlet_props, p_out, eta).unwrap_err();

        match err {
            CompressionError::ThermodynamicModelFailed { context, source } => {
                assert_eq!(context, "enthalpy(ideal outlet)");
                let source = source
                    .downcast_ref::<PropertyError>()
                    .expect("expected PropertyError source");
                assert!(matches!(source, PropertyError::Calculation { .. }));
            }
            other => panic!("unexpected error: {other:?}"),
        }
    }

    #[test]
    fn core_non_physical_work_returns_outlet_state() {
        let thermo = FakeThermo {
            mode: FakeMode::FixedEnthalpy(enth_si(-1.0)),
        };
        let (p_in, h_in, s_in, p_out, eta) = core_inputs();

        let inlet_props = InletProperties {
            thermo: &thermo,
            fluid: MockGas,
            p_in,
            h_in,
            s_in,
        };
        let err = isentropic_core(inlet_props, p_out, eta).unwrap_err();

        match err {
            CompressionError::NonPhysicalWork {
                outlet: _,
                raw_work,
                ..
            } => {
                assert!(raw_work.value < 0.0);
            }
            other => panic!("unexpected error: {other:?}"),
        }
    }
}