enok/src/motor/iso52016.rs

//! Dynamisk timesberegning av netto energibehov (beregningspunkt A) etter
//! NS-EN ISO 52016-1, slik NS 3031:2025 krever.
//!
//! Hver opake bygningsdel modelleres med fem noder der varmekapasiteten fordeles etter
//! masseplassering (ISO 52016-1 tabell B.13/B.14), og vinduer med to noder uten masse.
//! Sonen har én luftnode med varmekapasitet for luft og inventar, og én strålingsnode
//! («stjernenode») som fordeler langbølget strålingsutveksling mellom innvendige overflater.
//! Systemet løses implisitt (bakover-Euler) for hver time.
//!
//! Varme- og kjølebehovet bestemmes som i standarden: systemet løses først uten tilført
//! effekt og deretter med en prøveeffekt. Siden systemet er lineært, gir interpolasjon
//! mellom de to løsningene nøyaktig den effekten som treffer settpunktet.

use crate::klima::{sol, Vaerdata};
use crate::modell::lost::{Lost, Regulering};
use crate::modell::Mot;
use crate::motor::{lin, ventilasjon, RHO_C_LUFT};
use crate::{Result, TIMER};
use serde::Serialize;

/// Konvektiv varmeovergang innvendig overflate, W/(m²K). Gir sammen med h_ri R_si = 0,13.
const H_CI: f64 = 2.5;
/// Langbølget strålingsutveksling innvendig, W/(m²K).
const H_RI: f64 = 5.13;
/// Konvektiv varmeovergang utvendig, W/(m²K).
const H_CE: f64 = 20.0;
/// Langbølget stråling utvendig, W/(m²K). Gir sammen med h_ce R_se = 0,04.
const H_RE: f64 = 4.14;
/// Konvektiv andel av interne varmetilskudd (ISO 52016-1).
const KONVEKTIV_ANDEL: f64 = 0.4;
/// Innstråling der solskjerming trer i kraft, W/m².
const SKJERMING_TERSKEL: f64 = 300.0;
/// Utetemperatur der manuell solskjerming antas tatt i bruk, °C.
const SKJERMING_TEMP: f64 = 10.0;
const DT: f64 = 3600.0;
/// Prøveeffekt per m² ved bestemmelse av varme-/kjølebehov.
const PROVEEFFEKT: f64 = 10.0;

/// Nodeoppsett for én bygningsdel.
struct Del {
    areal: f64,
    /// Indeks til delens første node i tilstandsvektoren.
    forste: usize,
    /// Antall noder (5 for opake deler, 2 for vinduer og dører).
    antall: usize,
    /// Konduktans mellom nabonoder, W/K (allerede multiplisert med areal).
    h_indre: f64,
    /// Konduktans mot ytre miljø fra node 1, W/K.
    h_ytre: f64,
    /// Varmekapasitet per node, J/K.
    kapasitet: Vec<f64>,
    glass: bool,
    mot: Mot,
    sone: Option<usize>,
    b_faktor: f64,
    absorptans: f64,
    helning: f64,
    /// Solinnstråling på flaten, W/m², for hele året.
    innstraaling: Vec<f64>,
    g_effektiv: f64,
    skjermingsfaktor: f64,
    automatisk_skjerming: bool,
    horisontfaktor: f64,
}

/// Resultatet av en årssimulering: alle timesserier, W (energistrømmer) og °C.
#[derive(Serialize, Clone, Debug)]
pub struct Timeserie {
    pub t_ute: Vec<f64>,
    pub t_inne: Vec<f64>,
    pub t_operativ: Vec<f64>,
    /// Romoppvarming, post 1a, W.
    pub romoppvarming: Vec<f64>,
    /// Ventilasjonsvarme, post 1b, W.
    pub ventilasjonsvarme: Vec<f64>,
    /// Romkjøling, post 3a, W.
    pub romkjoling: Vec<f64>,
    /// Ventilasjonskjøling, post 3b, W.
    pub ventilasjonskjoling: Vec<f64>,
    /// Varmtvann, post 2, W.
    pub varmtvann: Vec<f64>,
    /// Vifter, post 4a, W.
    pub vifter: Vec<f64>,
    /// Pumper, post 4b, W.
    pub pumper: Vec<f64>,
    /// Belysning, post 5, W.
    pub belysning: Vec<f64>,
    /// Teknisk utstyr, post 6, W.
    pub utstyr: Vec<f64>,
    /// Sum solinnstråling gjennom vinduer, W.
    pub soltilskudd: Vec<f64>,
    /// Interne varmetilskudd, W.
    pub internlast: Vec<f64>,
    /// Solstråling på horisontal flate, W/m².
    pub sol_horisontal: Vec<f64>,
    /// Luftmengde ved vinduslufting, m³/h.
    pub vinduslufting: Vec<f64>,
    /// Temperatur i hver uoppvarmede sone, °C.
    pub sonetemperatur: Vec<Vec<f64>>,
}

impl Timeserie {
    fn ny() -> Self {
        let n = || vec![0.0; TIMER];
        Timeserie {
            t_ute: n(),
            t_inne: n(),
            t_operativ: n(),
            romoppvarming: n(),
            ventilasjonsvarme: n(),
            romkjoling: n(),
            ventilasjonskjoling: n(),
            varmtvann: n(),
            vifter: n(),
            pumper: n(),
            belysning: n(),
            utstyr: n(),
            soltilskudd: n(),
            internlast: n(),
            sol_horisontal: n(),
            vinduslufting: n(),
            sonetemperatur: Vec::new(),
        }
    }

    /// Årssum av en serie, kWh.
    pub fn sum_kwh(serie: &[f64]) -> f64 {
        serie.iter().sum::<f64>() / 1000.0
    }

    /// Månedssummer, kWh.
    pub fn maaned_kwh(serie: &[f64]) -> [f64; 12] {
        let mut ut = [0.0; 12];
        for (t, v) in serie.iter().enumerate() {
            ut[crate::maaned_for_time(t)] += v / 1000.0;
        }
        ut
    }
}

/// Uregulerbart varmetilskudd i sonen fra systemtap, W.
///
/// Tap fra tappevannsbereder og varmtvannsrør i oppvarmet sone varmer rommene og reduserer
/// romoppvarmingsbehovet, samtidig som energien telles under tappevann. Distribusjonstap for
/// romvarme i oppvarmet sone regnes ikke med her: den varmen ville uansett blitt avgitt i sonen,
/// og påvirker derfor bare fordelingen mellom rør og varmegiver, ikke energibehovet.
#[derive(Clone, Debug, Default)]
pub struct Systemtilskudd {
    /// Konstant tilskudd fra tappevannsanlegget i den oppvarmede sonen, W.
    pub tappevann: f64,
    /// Fast varmetilskudd i hver uoppvarmede sone, W.
    pub sone_fast: Vec<f64>,
    /// Varmetilskudd i hver uoppvarmede sone som andel av romoppvarmingseffekten.
    /// Dette er distribusjonstapet fra varmerør som går gjennom sonen.
    pub sone_andel_romvarme: Vec<f64>,
}

pub fn simuler(l: &Lost, v: &Vaerdata, tilskudd: Systemtilskudd) -> Result<Timeserie> {
    let deler = bygg_deler(l, v);
    let n_noder: usize = 2 + deler.iter().map(|d| d.antall).sum::<usize>();
    let areal_innvendig: f64 = deler.iter().map(|d| d.areal).sum();
    let c_luft = l.innvendig_varmekapasitet * l.bra;
    let h_inf = ventilasjon::infiltrasjonsluftmengde(l) * RHO_C_LUFT;
    let bakketemp = bakketemperatur(v);

    let mut ts = Timeserie::ny();
    ts.sonetemperatur = vec![vec![0.0; TIMER]; l.soner.len()];
    let mut t_soner: Vec<f64> = l
        .soner
        .iter()
        .map(|_| l.settpunkt_oppvarming - 5.0)
        .collect();
    let mut romvarme_forrige = 0.0;
    let mut tilstand = vec![l.settpunkt_oppvarming; n_noder];
    let mut t_a_forrige = l.settpunkt_oppvarming;

    // Innkjøring: gjenta de første to ukene til starttilstanden har stabilisert seg.
    for innkjoring in 0..3 {
        let skriv = innkjoring == 2;
        let start_tilstand = tilstand.clone();
        let timer: Box<dyn Iterator<Item = usize>> = if skriv {
            Box::new(0..TIMER)
        } else {
            Box::new(0..(24 * 14))
        };
        if !skriv {
            tilstand = start_tilstand;
        }
        for t in timer {
            let klokke = t % 24;
            let driftsdag = l.laster.driftsdag(t);
            let t_ute = v.t_ute[t];
            let t_himmel = sol::himmeltemperatur(v, t);

            // ── laster ──────────────────────────────────────────────
            // Lastprofilene gjelder driftsdager; på andre dager er bygningen ute av bruk.
            let skala = if driftsdag { l.bra } else { 0.0 };
            let bel = l.laster.belysning[klokke] * skala;
            let utstyr_el = l.laster.utstyr_effekt[klokke] * skala;
            let utstyr_varme = l.laster.utstyr_varme[klokke] * skala;
            let personer = l.laster.personer[klokke] * skala;
            let vv = l.laster.varmtvann[klokke] * skala;
            let pumper = l.laster.pumper[klokke] * skala;
            let vifter = ventilasjon::vifteeffekt(l, t);
            // Vifteeffekten ender som varme i lufta (ISO 52016-1: hele SFP-effekten).
            let internlast = bel + utstyr_varme + personer + tilskudd.tappevann + vifter;

            // ── sol ─────────────────────────────────────────────────
            let mut sol_gjennom = 0.0;
            for d in &deler {
                if d.glass && d.g_effektiv > 0.0 {
                    sol_gjennom += soltransmisjon(d, t, t_ute);
                }
            }

            // ── ventilasjon ─────────────────────────────────────────
            let q_vent = ventilasjon::luftmengde(l, t);
            let h_ve = q_vent * RHO_C_LUFT;
            let t_sup_hx = ventilasjon::tilluftstemperatur(l, t_ute, t_a_forrige);
            let kjol_tilluft = l
                .ventilasjon
                .tilluft_kjoling
                .filter(|&k| ventilasjon::i_drift(l, t) && t_sup_hx > k);
            let (t_sup, vent_varme, vent_kjol) =
                if l.ventilasjon.varmebatteri && t_sup_hx < l.ventilasjon.tilluft_settpunkt {
                    (
                        l.ventilasjon.tilluft_settpunkt,
                        h_ve * (l.ventilasjon.tilluft_settpunkt - t_sup_hx),
                        0.0,
                    )
                } else if let Some(k) = kjol_tilluft {
                    (k, 0.0, h_ve * (t_sup_hx - k))
                } else {
                    (t_sup_hx, 0.0, 0.0)
                };

            // ── settpunkt ───────────────────────────────────────────
            let sett_varme = if driftsdag && l.laster.driftstider.oppvarming.aktiv(klokke) {
                l.settpunkt_oppvarming
            } else {
                l.settpunkt_natt
            };

            // ── uoppvarmede soner ───────────────────────────────────
            // Temperaturen følger av den stasjonære varmebalansen for sonen, med
            // varmetilskudd fra rør og bereder som står der. Romoppvarmingen fra forrige
            // time brukes, siden rørtapet følger varmeanleggets drift.
            let t_bakke_naa = bakketemp[crate::maaned_for_time(t)];
            for (i, sone) in l.soner.iter().enumerate() {
                let tilskudd = tilskudd.sone_fast.get(i).copied().unwrap_or(0.0)
                    + tilskudd.sone_andel_romvarme.get(i).copied().unwrap_or(0.0)
                        * romvarme_forrige;
                t_soner[i] = sone.temperatur(t_a_forrige, t_ute, t_bakke_naa, tilskudd);
            }

            // ── løsning ─────────────────────────────────────────────
            let prove = PROVEEFFEKT * l.bra;
            let ctrl = |s: &[f64]| -> f64 {
                match l.regulering {
                    Regulering::Luft => s[0],
                    Regulering::Operativ => 0.5 * s[0] + 0.5 * s[1],
                }
            };
            let lag_ctx = |h_luft: f64| Kontekst {
                deler: &deler,
                n_noder,
                areal_innvendig,
                c_luft,
                h_inf,
                h_ve: h_ve + h_luft,
                // Lufteluften kommer rett utenfra, uten varmegjenvinning.
                t_sup: if h_ve + h_luft > 1e-9 {
                    (h_ve * t_sup + h_luft * t_ute) / (h_ve + h_luft)
                } else {
                    t_sup
                },
                t_ute,
                t_himmel,
                t_bakke: t_bakke_naa,
                t_a_forrige,
                t_soner: &t_soner,
                internlast,
                sol_gjennom,
                time: t,
            };

            // Vinduslufting: øk luftmengden til den operative temperaturen holder settpunktet.
            let mut h_luft = 0.0;
            let mut ctx = lag_ctx(0.0);
            let mut los0 = ctx.los(&tilstand, 0.0)?;
            if let Some((sett_luft, maks_luft)) = l.lufting {
                if ctrl(&los0) > sett_luft {
                    let h_maks = maks_luft * l.bra * RHO_C_LUFT;
                    let ctx_maks = lag_ctx(h_maks);
                    let los_maks = ctx_maks.los(&tilstand, 0.0)?;
                    if ctrl(&los_maks) >= sett_luft {
                        h_luft = h_maks;
                        ctx = ctx_maks;
                        los0 = los_maks;
                    } else {
                        let (mut lav, mut hoy) = (0.0, h_maks);
                        for _ in 0..12 {
                            let midt = 0.5 * (lav + hoy);
                            let k = lag_ctx(midt);
                            let l0 = k.los(&tilstand, 0.0)?;
                            if ctrl(&l0) > sett_luft {
                                lav = midt;
                            } else {
                                hoy = midt;
                            }
                            h_luft = midt;
                            ctx = k;
                            los0 = l0;
                        }
                    }
                }
            }
            let los1 = ctx.los(&tilstand, prove)?;
            let c0 = ctrl(&los0);
            let c1 = ctrl(&los1);
            let stigning = (c1 - c0) / prove; // K per W

            let mut fi = 0.0;
            if c0 < sett_varme - 1e-9 && stigning > 1e-12 {
                fi = (sett_varme - c0) / stigning;
                if let Some(maks) = l.maks_varmeeffekt {
                    fi = fi.min(maks);
                }
            } else if let Some(sett_kjol) = l.settpunkt_kjoling {
                if c0 > sett_kjol + 1e-9 && stigning > 1e-12 {
                    fi = (sett_kjol - c0) / stigning;
                    if let Some(maks) = l.maks_kjoleeffekt {
                        fi = fi.max(-maks);
                    }
                }
            }
            // Eksakt tilstand ved interpolasjon (systemet er lineært).
            let andel = fi / prove;
            for i in 0..n_noder {
                tilstand[i] = los0[i] + andel * (los1[i] - los0[i]);
            }
            t_a_forrige = tilstand[0];
            romvarme_forrige = fi.max(0.0);

            if skriv {
                ts.t_ute[t] = t_ute;
                ts.t_inne[t] = tilstand[0];
                ts.t_operativ[t] = 0.5 * tilstand[0] + 0.5 * tilstand[1];
                ts.romoppvarming[t] = fi.max(0.0);
                ts.romkjoling[t] = (-fi).max(0.0);
                ts.ventilasjonsvarme[t] = vent_varme;
                ts.ventilasjonskjoling[t] = vent_kjol;
                ts.varmtvann[t] = vv;
                ts.vifter[t] = vifter;
                ts.pumper[t] = pumper;
                ts.belysning[t] = bel;
                ts.utstyr[t] = utstyr_el;
                ts.soltilskudd[t] = sol_gjennom;
                ts.internlast[t] = internlast;
                ts.sol_horisontal[t] = v.ghi[t];
                ts.vinduslufting[t] = h_luft / RHO_C_LUFT;
                for (i, v) in t_soner.iter().enumerate() {
                    ts.sonetemperatur[i][t] = *v;
                }
            }
        }
    }
    Ok(ts)
}

struct Kontekst<'a> {
    deler: &'a [Del],
    n_noder: usize,
    areal_innvendig: f64,
    c_luft: f64,
    h_inf: f64,
    h_ve: f64,
    t_sup: f64,
    t_ute: f64,
    t_himmel: f64,
    t_bakke: f64,
    t_a_forrige: f64,
    t_soner: &'a [f64],
    internlast: f64,
    sol_gjennom: f64,
    time: usize,
}

impl Kontekst<'_> {
    /// Setter opp og løser timesligningene med gitt tilført effekt til sonen.
    fn los(&self, forrige: &[f64], fi_hc: f64) -> Result<Vec<f64>> {
        let n = self.n_noder;
        let mut a = vec![0.0; n * n];
        let mut b = vec![0.0; n];
        let idx = |i: usize, j: usize| i * n + j;

        // ── rad 0: sonens luftnode ──────────────────────────────────
        let c_dt = self.c_luft / DT;
        a[idx(0, 0)] = c_dt + self.h_ve + self.h_inf;
        b[0] = c_dt * forrige[0]
            + self.h_ve * self.t_sup
            + self.h_inf * self.t_ute
            + KONVEKTIV_ANDEL * self.internlast
            + fi_hc;
        for d in self.deler {
            let indre = d.forste + d.antall - 1;
            let h = d.areal * H_CI;
            a[idx(0, 0)] += h;
            a[idx(0, indre)] -= h;
        }

        // ── rad 1: strålingsnoden (ingen masse, netto null) ─────────
        for d in self.deler {
            let indre = d.forste + d.antall - 1;
            let h = d.areal * H_RI;
            a[idx(1, 1)] += h;
            a[idx(1, indre)] -= h;
        }
        b[1] = 0.0;

        // ── bygningsdelene ──────────────────────────────────────────
        let rad_last = (1.0 - KONVEKTIV_ANDEL) * self.internlast + self.sol_gjennom;
        for d in self.deler {
            let t_ytre = self.ytre_temperatur(d);
            let sol_absorbert = if d.glass {
                0.0
            } else {
                d.areal * d.absorptans * d.innstraaling[self.time]
            };
            for k in 0..d.antall {
                let i = d.forste + k;
                let c = d.kapasitet[k] / DT;
                a[idx(i, i)] += c;
                b[i] += c * forrige[i];
                if k == 0 {
                    a[idx(i, i)] += d.h_ytre;
                    b[i] += d.h_ytre * t_ytre + sol_absorbert;
                }
                if k + 1 < d.antall {
                    a[idx(i, i)] += d.h_indre;
                    a[idx(i, i + 1)] -= d.h_indre;
                }
                if k > 0 {
                    a[idx(i, i)] += d.h_indre;
                    a[idx(i, i - 1)] -= d.h_indre;
                }
                if k == d.antall - 1 {
                    a[idx(i, i)] += d.areal * (H_CI + H_RI);
                    a[idx(i, 0)] -= d.areal * H_CI;
                    a[idx(i, 1)] -= d.areal * H_RI;
                    b[i] += rad_last * d.areal / self.areal_innvendig;
                }
            }
        }
        lin::los(&mut a, &mut b, n)?;
        Ok(b)
    }

    fn ytre_temperatur(&self, d: &Del) -> f64 {
        match d.mot {
            Mot::Grunn => self.t_bakke,
            Mot::Uoppvarmet => match d.sone.and_then(|i| self.t_soner.get(i)) {
                Some(t) => *t,
                None => self.t_a_forrige - d.b_faktor * (self.t_a_forrige - self.t_ute),
            },
            Mot::Ute => {
                // Himmelandel etter flatens helning; resten ser bakken med utelufttemperatur.
                let f_himmel = (1.0 + d.helning.to_radians().cos()) / 2.0;
                let t_str = f_himmel * self.t_himmel + (1.0 - f_himmel) * self.t_ute;
                (H_CE * self.t_ute + H_RE * t_str) / (H_CE + H_RE)
            }
        }
    }
}

fn soltransmisjon(d: &Del, t: usize, t_ute: f64) -> f64 {
    let i_sol = d.innstraaling[t] * d.horisontfaktor;
    let skjerm = if d.skjermingsfaktor < 1.0
        && i_sol > SKJERMING_TERSKEL
        && (d.automatisk_skjerming || t_ute > SKJERMING_TEMP)
    {
        d.skjermingsfaktor
    } else {
        1.0
    };
    d.areal * d.g_effektiv * skjerm * i_sol
}

fn bygg_deler(l: &Lost, v: &Vaerdata) -> Vec<Del> {
    let mut deler = Vec::with_capacity(l.flater.len());
    let mut neste = 2usize;
    for f in &l.flater {
        let antall = if f.glass { 2 } else { 5 };
        let r_total = 1.0 / f.u_verdi.max(1e-6);
        let r_konstruksjon = (r_total - 1.0 / (H_CI + H_RI) - 1.0 / (H_CE + H_RE)).max(1e-3);
        let gap = (antall - 1) as f64;
        let h_indre = f.areal * gap / r_konstruksjon;
        let h_ytre = f.areal * (H_CE + H_RE);
        let kapasitet = fordel_kapasitet(f, antall);
        let innstraaling = if f.glass || f.mot == Mot::Ute {
            sol::innstraaling_aar(v, f.helning, f.orientering, l.albedo)
        } else {
            vec![0.0; TIMER]
        };
        deler.push(Del {
            areal: f.areal,
            forste: neste,
            antall,
            h_indre,
            h_ytre,
            kapasitet,
            glass: f.glass,
            mot: f.mot,
            sone: f.sone,
            b_faktor: f.b_faktor,
            absorptans: f.absorptans,
            helning: f.helning,
            innstraaling,
            g_effektiv: f.g_effektiv,
            skjermingsfaktor: f.skjermingsfaktor,
            automatisk_skjerming: f.automatisk_skjerming,
            horisontfaktor: (1.0 - f.horisont / 90.0).clamp(0.0, 1.0),
        });
        neste += antall;
    }
    deler
}

/// Fordeler varmekapasiteten over nodene etter ISO 52016-1 tabell B.13.
fn fordel_kapasitet(f: &crate::modell::lost::LostFlate, antall: usize) -> Vec<f64> {
    use crate::modell::Masseplassering::*;
    if antall == 2 {
        return vec![0.0, 0.0];
    }
    let k = f.kappa * f.areal;
    let mut c = vec![0.0; 5];
    match f.masseplassering {
        Innvendig => c[4] = k,
        Utvendig => c[0] = k,
        Begge => {
            c[0] = k / 2.0;
            c[4] = k / 2.0;
        }
        Midt => c[2] = k,
        Fordelt => {
            c[0] = k / 8.0;
            c[1] = k / 4.0;
            c[2] = k / 4.0;
            c[3] = k / 4.0;
            c[4] = k / 8.0;
        }
    }
    c
}

/// Dempet bakketemperatur per måned etter NS-EN ISO 13370 (forenklet): årsmiddel med
/// halvert amplitude og én måneds forsinkelse.
fn bakketemperatur(v: &Vaerdata) -> [f64; 12] {
    let aar = v.aarsmiddeltemperatur();
    let maaned = v.maanedsmiddeltemperatur();
    let mut ut = [0.0; 12];
    for m in 0..12 {
        let forsinket = maaned[(m + 11) % 12];
        ut[m] = aar + 0.5 * (forsinket - aar);
    }
    ut
}