simulatore-induttori/simulate.py

import csv
import math
import random
import shutil
from copy import deepcopy
from pathlib import Path

import numpy as np

from config import ARIA, PIASTRA, RANDOMIZZAZIONE, SENSORE, SIMULAZIONE, SORGENTE
from materials import MATERIALI


MU0 = 4.0 * math.pi * 1e-7


def calcola_skin_depth_m(materiale: dict, frequenza_hz: float) -> float:
    # Skin depth elettromagnetica approssimata:
    # delta = sqrt(2 * rho_e / (omega * mu))
    #
    # Semplificata. Per acciai ferromagnetici il comportamento reale
    # è fortemente non lineare con temperatura e campo magnetico.
    rho_e = materiale["resistivita_elettrica_ohm_m"]
    mu_r = materiale["permeabilita_relativa"]
    omega = 2.0 * math.pi * frequenza_hz
    mu = MU0 * mu_r
    return math.sqrt(2.0 * rho_e / (omega * mu))


def _spread_sorgenti_m(sorgente: dict) -> float:
    # Distanza lungo x tra la prima e l'ultima sorgente del gruppo.
    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
    return (numero_sorgenti - 1) * distanza


def _x_riferimento_iniziale_m(sorgente: dict) -> float:
    # Posizione a t=0 della sorgente di indice 0 (quella più arretrata nel
    # verso di marcia). x_inizio_m è la distanza dal sensore della sorgente
    # più avanzata (che quindi lo raggiunge per prima).
    spread = _spread_sorgenti_m(sorgente)
    x_sensore = sorgente["x_sensore_m"]
    x_inizio = sorgente["x_inizio_m"]
    if sorgente["velocita_m_s"] >= 0:
        return (x_sensore - x_inizio) - spread
    return (x_sensore + x_inizio) + spread


def _x_riferimento_finale_m(sorgente: dict) -> float:
    # Posizione di fine corsa della sorgente di indice 0. x_fine_m è la
    # distanza dal sensore della sorgente più arretrata (che quindi lo
    # supera per ultima).
    x_sensore = sorgente["x_sensore_m"]
    x_fine = sorgente["x_fine_m"]
    if sorgente["velocita_m_s"] >= 0:
        return x_sensore + x_fine
    return x_sensore - x_fine


def x_sorgente_al_tempo(sorgente: dict, t_s: float) -> float:
    return _x_riferimento_iniziale_m(sorgente) + sorgente["velocita_m_s"] * t_s


def _intervallo_attivo(inizio: float, fine: float, v: float, x_m: float) -> bool:
    if v >= 0:
        return inizio <= x_m <= fine
    return fine <= x_m <= inizio


def flusso_termico_incidente_W_m2(sorgente: dict, t_s: float) -> tuple[float, float]:
    # Restituisce x_sorgente_m (posizione della sorgente di riferimento) e
    # flusso_termico_efficace_W_m2 (somma dei contributi di tutte le sorgenti).
    #
    # Il movimento è rappresentato con un'impronta gaussiana centrata su
    # ciascuna sorgente in moto. Più sorgenti equidistanti si muovono insieme
    # come un gruppo rigido: condividono velocità, sigma e flusso di picco, e
    # sono sfalsate lungo x di un multiplo di "distanza_sorgenti_m". Il
    # modello 1D vede solo il flusso lungo la linea che passa per il sensore
    # fisso, sommando il contributo di tutte le sorgenti attive.
    x_rif_iniziale = _x_riferimento_iniziale_m(sorgente)
    x_rif_finale = _x_riferimento_finale_m(sorgente)
    x_riferimento = x_rif_iniziale + sorgente["velocita_m_s"] * t_s

    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
    v = sorgente["velocita_m_s"]
    zero_dopo_fine = sorgente.get("zero_dopo_fine", True)
    dy = sorgente["offset_y_percorso_m"]
    sigma = sorgente["sigma_punto_m"]

    q_totale = 0.0
    for i in range(numero_sorgenti):
        x_i = x_riferimento + i * distanza

        if zero_dopo_fine:
            inizio_i = x_rif_iniziale + i * distanza
            fine_i = x_rif_finale + i * distanza
            if not _intervallo_attivo(inizio_i, fine_i, v, x_i):
                continue

        dx = x_i - sorgente["x_sensore_m"]
        gaussiana = math.exp(-0.5 * (dx * dx + dy * dy) / (sigma * sigma))
        q_totale += (
            sorgente["flusso_termico_picco_W_m2"]
            * sorgente["efficienza_riscaldamento"]
            * gaussiana
        )

    return x_riferimento, q_totale


def riscaldamento_volumetrico_W_m3(
    q_superficie_W_m2: float,
    z_centri_m: np.ndarray,
    spessore_m: float,
    skin_depth_m: float,
) -> np.ndarray:
    # Converte il flusso superficiale equivalente in riscaldamento volumetrico q_vol(z).
    #
    # q_vol(z) = A * exp(-z / delta)
    #
    # A è scelto in modo che integrale_0^L q_vol(z) dz = q_superficie_W_m2.
    if q_superficie_W_m2 <= 0.0:
        return np.zeros_like(z_centri_m)

    delta = max(skin_depth_m, 1e-9)
    normalizzazione = delta * (1.0 - math.exp(-spessore_m / delta))
    return (q_superficie_W_m2 / normalizzazione) * np.exp(-z_centri_m / delta)


def costruisci_matrice_implicita(
    n: int,
    dt_s: float,
    dz_m: float,
    materiale: dict,
    h_caldo_W_m2K: float,
    h_freddo_W_m2K: float,
) -> np.ndarray:
    # Costruisce la matrice A per Eulero implicito:
    # A * T_next = rhs
    #
    # Le celle di bordo includono la convezione verso l'ambiente.
    k = materiale["conducibilita_termica_W_mK"]
    rho = materiale["densita_kg_m3"]
    cp = materiale["calore_specifico_J_kgK"]
    alpha = k / (rho * cp)

    r = alpha * dt_s / (dz_m * dz_m)
    b_caldo = h_caldo_W_m2K * dt_s / (rho * cp * dz_m)
    b_freddo = h_freddo_W_m2K * dt_s / (rho * cp * dz_m)

    A = np.zeros((n, n), dtype=float)

    # Bordo caldo, z = 0.
    A[0, 0] = 1.0 + r + b_caldo
    A[0, 1] = -r

    # Celle interne.
    for i in range(1, n - 1):
        A[i, i - 1] = -r
        A[i, i] = 1.0 + 2.0 * r
        A[i, i + 1] = -r

    # Bordo freddo, z = spessore.
    A[n - 1, n - 2] = -r
    A[n - 1, n - 1] = 1.0 + r + b_freddo

    return A


def quantizza(valore: float, passo: float) -> float:
    if passo <= 0.0:
        return valore
    return round(valore / passo) * passo


def configurazione_randomizzata(indice_run: int, rng: random.Random) -> dict:
    piastra = deepcopy(PIASTRA)
    aria = deepcopy(ARIA)
    sorgente = deepcopy(SORGENTE)
    sensore = deepcopy(SENSORE)

    if RANDOMIZZAZIONE.get("abilitata", False):
        def perturba_rel(valore: float, std_rel: float, fattore_min: float = 0.1) -> float:
            fattore = rng.gauss(1.0, std_rel)
            fattore = max(fattore_min, fattore)
            return valore * fattore

        sorgente["velocita_m_s"] = perturba_rel(
            sorgente["velocita_m_s"],
            RANDOMIZZAZIONE["velocita_std_rel"],
        )
        sorgente["flusso_termico_picco_W_m2"] = perturba_rel(
            sorgente["flusso_termico_picco_W_m2"],
            RANDOMIZZAZIONE["flusso_picco_std_rel"],
        )
        sorgente["sigma_punto_m"] = perturba_rel(
            sorgente["sigma_punto_m"],
            RANDOMIZZAZIONE["sigma_punto_std_rel"],
        )
        sorgente["offset_y_percorso_m"] = rng.uniform(
            -RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
            RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
        )
        aria["temperatura_ambiente_C"] += rng.gauss(
            0.0,
            RANDOMIZZAZIONE["temperatura_ambiente_std_C"],
        )
        sensore["rumore_std_C"] = perturba_rel(
            sensore["rumore_std_C"],
            RANDOMIZZAZIONE["rumore_sensore_std_rel"],
            fattore_min=0.0,
        )

    return {
        "id_run": f"run_{indice_run:04d}",
        "piastra": piastra,
        "aria": aria,
        "sorgente": sorgente,
        "sensore": sensore,
    }


def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
    piastra = cfg_run["piastra"]
    aria = cfg_run["aria"]
    sorgente = cfg_run["sorgente"]
    sensore = cfg_run["sensore"]

    nome_materiale = piastra["materiale"]
    materiale = MATERIALI[nome_materiale]

    spessore = piastra["spessore_mm"] / 1000.0
    n = piastra["n_nodi"]
    dz = spessore / n
    z_centri = (np.arange(n) + 0.5) * dz

    dt = SIMULAZIONE["dt_interno_s"]
    durata = SIMULAZIONE["durata_s"]
    periodo_campionamento = 1.0 / SIMULAZIONE["frequenza_campionamento_hz"]

    if sorgente["skin_depth_fissa_m"] is None:
        skin_depth = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
    else:
        skin_depth = float(sorgente["skin_depth_fissa_m"])

    A = costruisci_matrice_implicita(
        n=n,
        dt_s=dt,
        dz_m=dz,
        materiale=materiale,
        h_caldo_W_m2K=aria["h_caldo_W_m2K"],
        h_freddo_W_m2K=aria["h_freddo_W_m2K"],
    )

    rho = materiale["densita_kg_m3"]
    cp = materiale["calore_specifico_J_kgK"]

    b_caldo = aria["h_caldo_W_m2K"] * dt / (rho * cp * dz)
    b_freddo = aria["h_freddo_W_m2K"] * dt / (rho * cp * dz)

    T = np.full(n, piastra["temperatura_iniziale_C"], dtype=float)
    T_sensore = T[-1]

    prossimo_campione_t = 0.0
    T_vera_max = T[-1]
    T_misurata_max = T_sensore

    output_csv.parent.mkdir(parents=True, exist_ok=True)

    with output_csv.open("w", newline="") as f:
        writer = csv.writer(f)
        writer.writerow([
            "id_run",
            "tempo_s",
            "x_sorgente_m",
            "offset_y_sorgente_m",
            "flusso_termico_sorgente_W_m2",
            "skin_depth_m",
            "T_vera_lato_sensore_C",
            "T_misurata_sensore_C",
            "T_lato_caldo_C",
            "T_ambiente_C",
            "velocita_m_s",
            "sigma_punto_m",
            "flusso_picco_W_m2",
            "materiale",
        ])

        t = 0.0
        while t <= durata + 1e-12:
            x_sorgente, q_superficie = flusso_termico_incidente_W_m2(sorgente, t)
            q_vol = riscaldamento_volumetrico_W_m3(
                q_superficie_W_m2=q_superficie,
                z_centri_m=z_centri,
                spessore_m=spessore,
                skin_depth_m=skin_depth,
            )

            rhs = T + dt * q_vol / (rho * cp)
            rhs[0] += b_caldo * aria["temperatura_ambiente_C"]
            rhs[-1] += b_freddo * aria["temperatura_ambiente_C"]

            T = np.linalg.solve(A, rhs)

            # Temperatura vera sul lato freddo, dove si trova il sensore fisso.
            T_vera_lato_sensore = T[-1]

            # Inerzia del sensore del primo ordine.
            tau_sensore = max(sensore["costante_tempo_s"], 1e-9)
            T_sensore += (T_vera_lato_sensore - T_sensore) * dt / tau_sensore

            # Campionamento CSV.
            if t + 1e-12 >= prossimo_campione_t:
                misurata = T_sensore + rng.gauss(0.0, sensore["rumore_std_C"])
                misurata = quantizza(misurata, sensore["quantizzazione_C"])

                T_vera_max = max(T_vera_max, T_vera_lato_sensore)
                T_misurata_max = max(T_misurata_max, misurata)

                writer.writerow([
                    cfg_run["id_run"],
                    f"{t:.6f}",
                    f"{x_sorgente:.9f}",
                    f"{sorgente['offset_y_percorso_m']:.9f}",
                    f"{q_superficie:.6f}",
                    f"{skin_depth:.9e}",
                    f"{T_vera_lato_sensore:.6f}",
                    f"{misurata:.6f}",
                    f"{T[0]:.6f}",
                    f"{aria['temperatura_ambiente_C']:.6f}",
                    f"{sorgente['velocita_m_s']:.9f}",
                    f"{sorgente['sigma_punto_m']:.9f}",
                    f"{sorgente['flusso_termico_picco_W_m2']:.6f}",
                    nome_materiale,
                ])
                prossimo_campione_t += periodo_campionamento

            t += dt

    return {
        "id_run": cfg_run["id_run"],
        "file_csv": str(output_csv.name),
        "materiale": nome_materiale,
        "spessore_m": spessore,
        "n_nodi": n,
        "durata_s": durata,
        "frequenza_campionamento_hz": SIMULAZIONE["frequenza_campionamento_hz"],
        "dt_interno_s": dt,
        "temperatura_ambiente_C": aria["temperatura_ambiente_C"],
        "h_caldo_W_m2K": aria["h_caldo_W_m2K"],
        "h_freddo_W_m2K": aria["h_freddo_W_m2K"],
        "x_inizio_m": sorgente["x_inizio_m"],
        "x_fine_m": sorgente["x_fine_m"],
        "x_sensore_m": sorgente["x_sensore_m"],
        "offset_y_percorso_m": sorgente["offset_y_percorso_m"],
        "velocita_m_s": sorgente["velocita_m_s"],
        "numero_sorgenti": sorgente.get("numero_sorgenti", 1),
        "distanza_sorgenti_m": sorgente.get("distanza_sorgenti_m", 0.0),
        "sigma_punto_m": sorgente["sigma_punto_m"],
        "flusso_termico_picco_W_m2": sorgente["flusso_termico_picco_W_m2"],
        "efficienza_riscaldamento": sorgente["efficienza_riscaldamento"],
        "frequenza_hz": sorgente["frequenza_hz"],
        "skin_depth_m": skin_depth,
        "costante_tempo_sensore_s": sensore["costante_tempo_s"],
        "rumore_std_sensore_C": sensore["rumore_std_C"],
        "quantizzazione_sensore_C": sensore["quantizzazione_C"],
        "T_vera_max_lato_sensore_C": T_vera_max,
        "T_misurata_max_sensore_C": T_misurata_max,
    }


def main() -> None:
    cartella_output = Path(SIMULAZIONE["cartella_output"])
    if cartella_output.exists():
        shutil.rmtree(cartella_output)
    cartella_output.mkdir(parents=True, exist_ok=True)

    rng = random.Random(SIMULAZIONE["seed"])

    righe_metadata = []
    for i in range(1, SIMULAZIONE["num_run"] + 1):
        cfg_run = configurazione_randomizzata(i, rng)
        percorso_csv = cartella_output / f"{cfg_run['id_run']}.csv"
        righe_metadata.append(simula_singolo(cfg_run, percorso_csv, rng))

    percorso_metadata = cartella_output / "metadata.csv"
    with percorso_metadata.open("w", newline="") as f:
        writer = csv.DictWriter(f, fieldnames=list(righe_metadata[0].keys()))
        writer.writeheader()
        writer.writerows(righe_metadata)

    print(f"Generati {len(righe_metadata)} run in: {cartella_output.resolve()}")
    print(f"Metadata: {percorso_metadata.resolve()}")


if __name__ == "__main__":
    main()