simulatore-induttori/fem.py

# Solutore termico transitorio a elementi finiti sulla mesh shell della
# fascetta (vedi mesh.py per la geometria).
#
# Il dominio è la superficie media del cilindro, con x e la coordinata
# circonferenziale s = R·theta entrambe risolte, mentre lo spessore è
# collassato (temperatura uniforme attraverso la parete).
#
# Il collasso in spessore è lecito: con spessore 0.18 mm il numero di Biot
# h·t/k vale ~1e-7 e il tempo di diffusione attraverso la parete t²/alpha è
# di pochi millisecondi, molto più rapido del transito delle sorgenti. Di
# conseguenza la profondità di penetrazione del riscaldamento a induzione non
# entra nel modello: conta solo il flusso totale assorbito per unità di
# superficie. La skin depth è comunque calcolata dalla frequenza dell'induttore
# come grandezza diagnostica (`calcola_skin_depth_m`).
#
# Equazione risolta (per unità di superficie media):
#
#   rho·cp·t ∂T/∂t = k·t ∇²T + q(x, s, t) - (h_est + h_int)·(T - T_amb)
#
# con convezione aggiuntiva sui due bordi anulari x = 0 e x = L (coefficiente
# h_bordi su un'area pari a spessore × perimetro). L'impronta della sorgente
# è una gaussiana isotropa nel piano (x, s), con la distanza circonferenziale
# valutata sull'immagine più vicina perché la superficie è chiusa.
#
# Discretizzazione: elementi shell quadrangolari bilineari a 4 nodi. Gli
# elementi della mesh sono tutti rettangoli identici dx × ds, quindi le
# matrici di elemento sono calcolate una volta sola e assemblate in forma
# vettorizzata. Integrazione temporale con Eulero implicito: la matrice di
# sistema è costante e viene fattorizzata LU una volta per run.

import math
import random
from copy import deepcopy

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import splu

from config import ARIA, FASCETTA, FEM, RANDOMIZZAZIONE, SENSORE, SORGENTE
from materials import MATERIALI
from mesh import genera_mesh

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

# Matrici di riferimento dell'elemento rettangolare bilineare a 4 nodi, con
# nodi locali in ordine antiorario: (0,0), (a,0), (a,b), (0,b).
#
# MASSA_RIF va moltiplicata per l'area a·b e dà l'integrale di N_i·N_j;
# RIGIDEZZA_X per b/a e RIGIDEZZA_S per a/b, e la loro somma dà l'integrale
# di grad(N_i)·grad(N_j).
MASSA_RIF = np.array(
    [[4.0, 2.0, 1.0, 2.0],
     [2.0, 4.0, 2.0, 1.0],
     [1.0, 2.0, 4.0, 2.0],
     [2.0, 1.0, 2.0, 4.0]]
) / 36.0

RIGIDEZZA_X = np.array(
    [[2.0, -2.0, -1.0, 1.0],
     [-2.0, 2.0, 1.0, -1.0],
     [-1.0, 1.0, 2.0, -2.0],
     [1.0, -1.0, -2.0, 2.0]]
) / 6.0

RIGIDEZZA_S = np.array(
    [[2.0, 1.0, -1.0, -2.0],
     [1.0, 2.0, -2.0, -1.0],
     [-1.0, -2.0, 2.0, 1.0],
     [-2.0, -1.0, 1.0, 2.0]]
) / 6.0


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.
    #
    # Nel modello shell è una grandezza diagnostica: serve a verificare che la
    # deposizione del calore resti confinata entro uno spessore confrontabile
    # con la parete, non entra nell'equazione risolta.
    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 quantizza(valore: float, passo: float) -> float:
    if passo <= 0.0:
        return valore
    return round(valore / passo) * passo


def configurazione_randomizzata(rng: random.Random | None = None) -> dict:
    """Copia dei dizionari di configurazione usati da un'analisi FEM.

    Se `RANDOMIZZAZIONE["abilitata"]`, velocità, flusso di picco, sigma dello
    spot, offset y del percorso, temperatura ambiente e rumore del sensore sono
    perturbati con estrazioni dall'RNG passato: analisi diverse dello stesso
    setup differiscono come differirebbero due passaggi reali.
    """
    fascetta = deepcopy(FASCETTA)
    aria = deepcopy(ARIA)
    sorgente = deepcopy(SORGENTE)
    sensore = deepcopy(SENSORE)

    if RANDOMIZZAZIONE.get("abilitata", False):
        if rng is None:
            rng = random.Random(FEM["seed"])

        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 {
        "fascetta": fascetta,
        "aria": aria,
        "sorgente": sorgente,
        "sensore": sensore,
    }


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, x_sensore_m: float) -> float:
    # Posizione a t=0 della sorgente di indice 0. Le sorgenti i sono a
    # x_i = riferimento + i * distanza, quindi per v >= 0 l'indice 0 è la
    # più arretrata nel verso di marcia, per v < 0 è la più avanzata.
    # x_inizio_m è la distanza dal sensore della sorgente più avanzata
    # (quella che lo raggiunge per prima).
    spread = _spread_sorgenti_m(sorgente)
    x_inizio = sorgente["x_inizio_m"]
    if sorgente["velocita_m_s"] >= 0:
        return (x_sensore_m - x_inizio) - spread
    return x_sensore_m + x_inizio


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


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 _assembla(elementi: np.ndarray, matrice_elemento: np.ndarray, n_nodi: int):
    # Assembla una matrice globale sparsa a partire da un'unica matrice di
    # elemento 4×4, uguale per tutti gli elementi. I contributi ripetuti sullo
    # stesso (riga, colonna) vengono sommati dal formato COO.
    n_elementi = len(elementi)
    righe = np.broadcast_to(elementi[:, :, None], (n_elementi, 4, 4)).ravel()
    colonne = np.broadcast_to(elementi[:, None, :], (n_elementi, 4, 4)).ravel()
    valori = np.broadcast_to(matrice_elemento, (n_elementi, 4, 4)).ravel()
    return sp.coo_matrix(
        (valori, (righe, colonne)), shape=(n_nodi, n_nodi)
    ).tocsr()


def _massa_anello(indici: np.ndarray, lunghezza_segmento_m: float, n_nodi: int):
    # Matrice di massa 1D su un anello chiuso di nodi (integrale di N_i·N_j
    # lungo la circonferenza), usata per la convezione sui bordi x = 0 e x = L.
    n = len(indici)
    successivi = indici[(np.arange(n) + 1) % n]
    coefficiente = lunghezza_segmento_m / 6.0
    righe = np.concatenate([indici, indici, successivi, successivi])
    colonne = np.concatenate([indici, successivi, indici, successivi])
    valori = coefficiente * np.concatenate(
        [np.full(n, 2.0), np.full(n, 1.0), np.full(n, 1.0), np.full(n, 2.0)]
    )
    return sp.coo_matrix(
        (valori, (righe, colonne)), shape=(n_nodi, n_nodi)
    ).tocsr()


def prepara_stato_fem(
    fascetta: dict,
    aria: dict,
    sorgente: dict,
    sensore: dict,
    mesh_dati: dict | None = None,
    dt_s: float | None = None,
) -> dict:
    """Assembla le matrici FEM e fattorizza il sistema implicito.

    Restituisce lo stato costante per l'integrazione temporale: matrici,
    solutore LU, coordinate nodali, peso circonferenziale della sorgente e
    indice del nodo osservato dal sensore.
    """
    if mesh_dati is None:
        mesh_dati = genera_mesh(fascetta)
    if dt_s is None:
        dt_s = FEM["dt_s"]

    materiale = MATERIALI[fascetta["materiale"]]
    if sorgente["skin_depth_fissa_m"] is None:
        skin_depth_m = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
    else:
        skin_depth_m = float(sorgente["skin_depth_fissa_m"])
    k = materiale["conducibilita_termica_W_mK"]
    rho = materiale["densita_kg_m3"]
    cp = materiale["calore_specifico_J_kgK"]

    spessore_m = mesh_dati["spessore_m"]
    raggio_m = mesh_dati["raggio_m"]
    circonferenza_m = 2.0 * math.pi * raggio_m
    n_x = mesh_dati["n_elementi_x"]
    n_theta = mesh_dati["n_elementi_circonferenza"]
    elementi = mesh_dati["elementi"]
    indice_nodo = mesh_dati["indice_nodo"]
    n_nodi = len(mesh_dati["nodi"])

    dx_m = mesh_dati["lunghezza_m"] / n_x
    ds_m = circonferenza_m / n_theta

    h_esterno = aria["h_esterno_W_m2K"]
    h_interno = aria["h_interno_W_m2K"]
    h_bordi = aria["h_bordi_W_m2K"]
    T_ambiente = aria["temperatura_ambiente_C"]

    # Matrice di massa di superficie: riusata per la capacità termica, per la
    # convezione sulle due facce e per il carico della sorgente.
    M_area = _assembla(elementi, MASSA_RIF * (dx_m * ds_m), n_nodi)
    K_cond = _assembla(
        elementi,
        k * spessore_m * (RIGIDEZZA_X * (ds_m / dx_m) + RIGIDEZZA_S * (dx_m / ds_m)),
        n_nodi,
    )

    # Convezione sui bordi anulari: l'area di scambio di ogni segmento è
    # spessore × lunghezza del segmento circonferenziale.
    H_bordi = h_bordi * spessore_m * (
        _massa_anello(indice_nodo[0, :], ds_m, n_nodi)
        + _massa_anello(indice_nodo[-1, :], ds_m, n_nodi)
    )

    C_su_dt = (rho * cp * spessore_m / dt_s) * M_area
    A = C_su_dt + K_cond + (h_esterno + h_interno) * M_area + H_bordi

    # Termine noto costante della convezione: (somma delle righe) × h × T_amb,
    # perché la somma delle righe della matrice di massa è l'integrale di N_i.
    area_nodale = np.asarray(M_area.sum(axis=1)).ravel()
    bordo_nodale = np.asarray(H_bordi.sum(axis=1)).ravel()
    carico_convezione = (h_esterno + h_interno) * area_nodale * T_ambiente
    carico_convezione += bordo_nodale * T_ambiente

    # Peso circonferenziale dell'impronta della sorgente: dipende solo
    # dall'offset y del percorso e dalla sigma dello spot, entrambi costanti
    # durante il run, quindi è valutato una volta sola qui. La distanza è
    # quella dell'immagine più vicina, perché la superficie è chiusa.
    sigma_m = sorgente["sigma_punto_m"]
    ds_sorgente = raggio_m * mesh_dati["theta_nodi_rad"] - sorgente["offset_y_percorso_m"]
    ds_sorgente -= circonferenza_m * np.round(ds_sorgente / circonferenza_m)
    peso_circonferenziale = np.exp(-0.5 * (ds_sorgente ** 2) / (sigma_m * sigma_m))

    # Nodo osservato dal sensore: x più vicino a x_mm sul piano theta = 0.
    x_sensore_m = sensore["x_mm"] / 1000.0
    i_sensore = int(np.argmin(np.abs(mesh_dati["x_nodi_m"] - x_sensore_m)))
    indice_sensore = int(indice_nodo[i_sensore, 0])

    return {
        "mesh": mesh_dati,
        "dt_s": dt_s,
        "M_area": M_area,
        "C_su_dt": C_su_dt,
        "carico_convezione": carico_convezione,
        "solutore": splu(sp.csc_matrix(A)),
        "x_nodi_m": mesh_dati["x_nodi_m"],
        "peso_circonferenziale": peso_circonferenziale,
        "forma_griglia": (n_x + 1, n_theta),
        "circonferenza_m": circonferenza_m,
        "area_nodale_m2": area_nodale,
        "T_ambiente_C": T_ambiente,
        "x_sensore_m": x_sensore_m,
        "indice_sensore": indice_sensore,
        "n_nodi": n_nodi,
        # Diagnostica: non entra nell'equazione risolta.
        "skin_depth_m": skin_depth_m,
    }


def flusso_nodale_W_m2(
    sorgente: dict,
    stato: dict,
    t_s: float,
) -> tuple[float, np.ndarray]:
    """Flusso termico assorbito nei nodi all'istante t_s.

    Restituisce (x_riferimento_m, q_nodale) dove x_riferimento_m è la
    posizione della sorgente di indice 0. Ogni sorgente attiva contribuisce
    con un'impronta gaussiana isotropa in (x, s).

    L'impronta gaussiana isotropa è separabile, q(x, s) = q_x(x) · q_s(s), e
    il fattore circonferenziale q_s è costante nel tempo (precalcolato in
    `prepara_stato_fem`). Qui si valuta quindi solo il fattore assiale sui
    n_x + 1 nodi distinti in x e si espande con un prodotto esterno, invece di
    valutare l'esponenziale su tutti i nodi della mesh.
    """
    x_sensore_m = stato["x_sensore_m"]
    x_rif_iniziale = _x_riferimento_iniziale_m(sorgente, x_sensore_m)
    x_rif_finale = _x_riferimento_finale_m(sorgente, x_sensore_m)
    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)
    sigma = sorgente["sigma_punto_m"]
    q_picco = sorgente["flusso_termico_picco_W_m2"] * sorgente["efficienza_riscaldamento"]

    peso_assiale = np.zeros(len(stato["x_nodi_m"]))
    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 = stato["x_nodi_m"] - x_i
        peso_assiale += np.exp(-0.5 * (dx * dx) / (sigma * sigma))

    q_nodale = q_picco * np.outer(peso_assiale, stato["peso_circonferenziale"])
    return x_riferimento, q_nodale.ravel()


def passo_implicito_fem(stato: dict, T: np.ndarray, q_nodale: np.ndarray) -> np.ndarray:
    # Avanza il campo nodale di un passo dt: termine noto = capacità sul campo
    # precedente + carico della sorgente + carico di convezione, poi risoluzione
    # del sistema già fattorizzato.
    rhs = stato["C_su_dt"] @ T
    rhs += stato["M_area"] @ q_nodale
    rhs += stato["carico_convezione"]
    return stato["solutore"].solve(rhs)


def campo_iniziale(stato: dict) -> np.ndarray:
    return np.full(stato["n_nodi"], stato["T_ambiente_C"], dtype=float)


def simula_campo_fem(
    dt_frame_s: float,
    cfg: dict | None = None,
    durata_s: float | None = None,
    mesh_dati: dict | None = None,
    rng: random.Random | None = None,
) -> dict:
    """Integra il campo FEM fino a durata_s, salvando un campo ogni dt_frame_s.

    Il campo nodale è la temperatura vera della parete; la lettura del sensore
    ne è la versione degradata: inerzia del primo ordine, rumore gaussiano e
    quantizzazione, questi ultimi due applicati a ogni campionamento.
    """
    if rng is None:
        rng = random.Random(FEM["seed"])
    if cfg is None:
        cfg = configurazione_randomizzata(rng)
    if durata_s is None:
        durata_s = FEM["durata_s"]

    fascetta = cfg["fascetta"]
    aria = cfg["aria"]
    sorgente = cfg["sorgente"]
    sensore = cfg["sensore"]

    stato = prepara_stato_fem(fascetta, aria, sorgente, sensore, mesh_dati=mesh_dati)
    dt = stato["dt_s"]
    indice_sensore = stato["indice_sensore"]

    T = campo_iniziale(stato)
    T_sensore = float(T[indice_sensore])
    tau_sensore = max(sensore["costante_tempo_s"], 1e-9)

    tempi, campi, x_riferimenti = [], [], []
    T_vere, T_lette = [], []

    prossimo_frame_t = 0.0
    t = 0.0
    while t <= durata_s + 1e-12:
        x_rif, q_nodale = flusso_nodale_W_m2(sorgente, stato, t)
        T = passo_implicito_fem(stato, T, q_nodale)

        T_sensore += (T[indice_sensore] - T_sensore) * dt / tau_sensore

        if t + 1e-12 >= prossimo_frame_t:
            letta = T_sensore + rng.gauss(0.0, sensore["rumore_std_C"])
            letta = quantizza(letta, sensore["quantizzazione_C"])

            tempi.append(t)
            campi.append(T.copy())
            x_riferimenti.append(x_rif)
            T_vere.append(float(T[indice_sensore]))
            T_lette.append(letta)
            prossimo_frame_t += dt_frame_s

        t += dt

    return {
        "stato": stato,
        "mesh": stato["mesh"],
        "tempi": np.array(tempi),
        "campi": campi,
        "x_riferimenti": np.array(x_riferimenti),
        "T_vere": np.array(T_vere),
        "T_lette": np.array(T_lette),
        "sorgente": sorgente,
        "T_ambiente_C": stato["T_ambiente_C"],
        "skin_depth_m": stato["skin_depth_m"],
    }