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 x_sorgente_al_tempo(sorgente: dict, t_s: float) -> float: return sorgente["x_inizio_m"] + sorgente["velocita_m_s"] * t_s def sorgente_attiva(sorgente: dict, x_m: float) -> bool: if not sorgente.get("zero_dopo_fine", True): return True inizio = sorgente["x_inizio_m"] fine = sorgente["x_fine_m"] v = sorgente["velocita_m_s"] 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 e flusso_termico_efficace_W_m2. # # Il movimento è rappresentato con un'impronta gaussiana centrata sulla # sorgente in moto. Il modello 1D vede solo il flusso lungo la linea # che passa per il sensore fisso. x = x_sorgente_al_tempo(sorgente, t_s) if not sorgente_attiva(sorgente, x): return x, 0.0 dx = x - sorgente["x_sensore_m"] dy = sorgente["offset_y_percorso_m"] sigma = sorgente["sigma_punto_m"] gaussiana = math.exp(-0.5 * (dx * dx + dy * dy) / (sigma * sigma)) q = ( sorgente["flusso_termico_picco_W_m2"] * sorgente["efficienza_riscaldamento"] * gaussiana ) return x, q 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_m"] 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"], "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()