From c67f6f8da63b8d57de399ee6788606268a1d3a45 Mon Sep 17 00:00:00 2001 From: Davide Grilli Date: Sun, 5 Jul 2026 19:22:51 +0200 Subject: [PATCH] Scambia calore su tutto il contorno e per conduzione circonferenziale - Convezione anche sui bordi laterali x = 0 e x = lunghezza (nuovo coefficiente h_bordi_W_m2K in ARIA) - Termine di conduzione lungo y verso il resto della fascetta, assunto a temperatura ambiente: q_y = -k*(T - T_amb)/L_y^2 con L_y configurabile (lunghezza_conduzione_y_mm in FASCETTA) - La temperatura iniziale del campo e' la temperatura ambiente del run (rimosso temperatura_iniziale_C) - Fisica estratta in prepara_stato_termico e passo_implicito, condivise tra simulate.py e plot_animazione.py per evitare duplicazione - Entrambi i termini sono lineari e restano nella matrice fattorizzata una volta per run Co-Authored-By: Claude Fable 5 --- CLAUDE.md | 9 ++- config.py | 13 +++- plot_animazione.py | 56 ++++------------- simulate.py | 153 ++++++++++++++++++++++++++++++--------------- 4 files changed, 130 insertions(+), 101 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index ea77e5f..f07e41f 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -47,10 +47,13 @@ Generatore di dataset per misurazioni termiche pseudo-realistiche di una fascett - La skin depth è calcolata dalla resistività elettrica del materiale e dalla frequenza di induzione (`calcola_skin_depth_m`) - Le sorgenti gaussiane in movimento producono un profilo di flusso termico superficiale q(x) sul lato esterno, variabile nel tempo (`profilo_flusso_incidente_W_m2`) - Quel flusso è ridistribuito volumetricamente attraverso lo spessore con decadimento esponenziale in z (`profilo_deposizione_z_1_m`): q_vol(x, z) = q(x) · p(z) -- Uno schema 2D a volumi finiti con Eulero implicito integra l'equazione del calore su `n_nodi_x × n_nodi_z` celle: la matrice sparsa è fattorizzata LU una volta per run (`costruisci_solutore_implicito_2d`, che restituisce l'oggetto `splu`), poi ad ogni passo si risolve solo il sistema triangolare -- Le condizioni al contorno sono incorporate nella matrice: convezione sui lati esterno (z = 0) e interno (z = spessore), bordi adiabatici a x = 0 e x = lunghezza +- Uno schema 2D a volumi finiti con Eulero implicito integra l'equazione del calore su `n_nodi_x × n_nodi_z` celle: `prepara_stato_termico` costruisce griglia, coefficienti e matrice sparsa fattorizzata LU una volta per run (`costruisci_solutore_implicito_2d`, che restituisce l'oggetto `splu`), poi `passo_implicito` avanza il campo risolvendo solo il sistema triangolare +- Le condizioni al contorno sono incorporate nella matrice: convezione su tutti e quattro i lati della sezione, più un termine di conduzione circonferenziale (y) verso il resto della fascetta assunto a temperatura ambiente (q_y = -k·(T - T_amb)/L_y², con L_y = `lunghezza_conduzione_y_mm`) +- La temperatura iniziale del campo è la temperatura ambiente (randomizzata per run) - L'output del sensore aggiunge inerzia del primo ordine, rumore gaussiano e quantizzazione +`prepara_stato_termico` e `passo_implicito` sono condivisi con `plot_animazione.py`, che riproduce la fisica di run_0001 per animare la sezione: ogni modifica alla fisica va fatta lì, non duplicata. + **Randomizzazione per run** (`configurazione_randomizzata`): ogni run perturba velocità, flusso di picco, sigma del punto, offset y, temperatura ambiente e rumore del sensore con estrazioni gaussiane/uniformi da un RNG con seed fisso, garantendo riproducibilità. **Schema di output** (`dataset/run_XXXX.csv`): serie temporale con colonne `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`. `metadata.csv` ha una riga per run con tutti i parametri e le temperature di picco. @@ -61,7 +64,7 @@ Ogni parametro in [config.py](config.py) ha un commento che spiega solo cos'è ( ## Vincoli progettuali chiave -- Il modello è 2D nella sezione (x = lunghezza, z = spessore). La coordinata circonferenziale y non è risolta spazialmente — l'offset y del percorso delle sorgenti è collassato in un'attenuazione gaussiana del flusso; il diametro è registrato solo come geometria del setup. +- Il modello è 2D nella sezione (x = lunghezza, z = spessore). La coordinata circonferenziale y non è risolta spazialmente — l'offset y del percorso delle sorgenti è collassato in un'attenuazione gaussiana del flusso, e la conduzione lungo y è un termine di scambio lineare verso la temperatura ambiente; il diametro è registrato solo come geometria del setup. - Le posizioni di inizio/fine corsa delle sorgenti (`x_inizio_m`, `x_fine_m`) sono distanze dal punto x del sensore lungo il verso di marcia; il segno di `velocita_m_s` determina il verso (negativo = -x). - La matrice implicita è costruita e fattorizzata una volta per run (proprietà del materiale costanti, nessun coefficiente dipendente dalla temperatura). Se si aggiungono proprietà dipendenti dalla temperatura, la matrice deve essere ricostruita e rifattorizzata ad ogni passo temporale. - `simulate.py` cancella e ricrea l'intera cartella di output ad ogni esecuzione (`shutil.rmtree`). diff --git a/config.py b/config.py index 88c0fb4..5eececb 100644 --- a/config.py +++ b/config.py @@ -14,6 +14,10 @@ # - Le sorgenti a induzione si muovono in direzione -x sul lato esterno. # - Il riscaldamento a induzione è approssimato come riscaldamento volumetrico # che decade esponenzialmente con la profondità z secondo la skin depth. +# - La sezione scambia per convezione con l'aria su tutti e quattro i lati +# (esterno, interno e i due bordi in x) e per conduzione lungo y con il +# resto della fascetta, assunto a temperatura ambiente. +# - La temperatura iniziale della fascetta è pari alla temperatura ambiente. # - Il sensore è un pirometro a infrarossi posto all'interno della fascetta, # a distanza "distanza_parete_mm" dalla parete interna: misura senza # contatto la temperatura della superficie interna nel punto x = "x_mm". @@ -64,8 +68,10 @@ FASCETTA = { # Più nodi = maggiore risoluzione spaziale, simulazione più lenta. "n_nodi_z": 15, - # Temperatura iniziale uniforme. - "temperatura_iniziale_C": 25.0, + # Distanza caratteristica lungo y (direzione circonferenziale) su cui la + # sezione scambia calore per conduzione con il resto della fascetta, + # assunto a temperatura ambiente. + "lunghezza_conduzione_y_mm": 25.0, # Deve corrispondere a una chiave in materials.py. "materiale": "banda_stagnata", @@ -80,6 +86,9 @@ ARIA = { # Coefficiente di convezione sul lato interno (z = spessore, lato sensore). "h_interno_W_m2K": 8.0, + + # Coefficiente di convezione sui bordi laterali (x = 0 e x = lunghezza). + "h_bordi_W_m2K": 10.0, } SORGENTE = { diff --git a/plot_animazione.py b/plot_animazione.py index c1d3d2a..be58e21 100644 --- a/plot_animazione.py +++ b/plot_animazione.py @@ -14,12 +14,10 @@ import numpy as np from matplotlib.animation import FuncAnimation, PillowWriter from config import SIMULAZIONE -from materials import MATERIALI from simulate import ( - calcola_skin_depth_m, configurazione_randomizzata, - costruisci_solutore_implicito_2d, - profilo_deposizione_z_1_m, + passo_implicito, + prepara_stato_termico, profilo_flusso_incidente_W_m2, ) @@ -43,45 +41,17 @@ def simula_campi(cfg_run: dict) -> dict: aria = cfg_run["aria"] sorgente = cfg_run["sorgente"] sensore = cfg_run["sensore"] - materiale = MATERIALI[fascetta["materiale"]] - lunghezza = fascetta["lunghezza_mm"] / 1000.0 - spessore = fascetta["spessore_mm"] / 1000.0 - n_x = fascetta["n_nodi_x"] - n_z = fascetta["n_nodi_z"] - dx = lunghezza / n_x - dz = spessore / n_z - x_centri = (np.arange(n_x) + 0.5) * dx - z_centri = (np.arange(n_z) + 0.5) * dz + stato = prepara_stato_termico(fascetta, aria, sorgente) + n_x = stato["n_x"] + n_z = stato["n_z"] + x_centri = stato["x_centri_m"] + dt = stato["dt_s"] x_sensore = sensore["x_mm"] / 1000.0 - i_sensore = min(n_x - 1, max(0, int(x_sensore / dx))) + i_sensore = min(n_x - 1, max(0, int(x_sensore / stato["dx_m"]))) - dt = SIMULAZIONE["dt_interno_s"] - - 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"]) - - solutore = costruisci_solutore_implicito_2d( - n_x=n_x, - n_z=n_z, - dt_s=dt, - dx_m=dx, - dz_m=dz, - materiale=materiale, - h_esterno_W_m2K=aria["h_esterno_W_m2K"], - h_interno_W_m2K=aria["h_interno_W_m2K"], - ) - - rho = materiale["densita_kg_m3"] - cp = materiale["calore_specifico_J_kgK"] - b_esterno = aria["h_esterno_W_m2K"] * dt / (rho * cp * dz) - b_interno = aria["h_interno_W_m2K"] * dt / (rho * cp * dz) - profilo_z = profilo_deposizione_z_1_m(z_centri, spessore, skin_depth) - - T = np.full((n_x, n_z), fascetta["temperatura_iniziale_C"], dtype=float) + T = np.full((n_x, n_z), aria["temperatura_ambiente_C"], dtype=float) T_sensore = T[i_sensore, -1] tau_sensore = max(sensore["costante_tempo_s"], 1e-9) @@ -92,11 +62,7 @@ def simula_campi(cfg_run: dict) -> dict: t = 0.0 while t <= T_FINE_ANIMAZIONE_S + 1e-12: x_rif, q_x = profilo_flusso_incidente_W_m2(sorgente, x_sensore, t, x_centri) - - rhs = T + (dt / (rho * cp)) * q_x[:, None] * profilo_z[None, :] - rhs[:, 0] += b_esterno * aria["temperatura_ambiente_C"] - rhs[:, -1] += b_interno * aria["temperatura_ambiente_C"] - T = solutore.solve(rhs.ravel()).reshape(n_x, n_z) + T = passo_implicito(stato, T, q_x) T_sensore += (T[i_sensore, -1] - T_sensore) * dt / tau_sensore @@ -176,7 +142,7 @@ def main() -> None: extent=(0.0, lunghezza_mm, spessore_mm, 0.0), aspect="auto", cmap="inferno", - vmin=cfg_run["fascetta"]["temperatura_iniziale_C"], + vmin=cfg_run["aria"]["temperatura_ambiente_C"], vmax=T_max, interpolation="bilinear", ) diff --git a/simulate.py b/simulate.py index 54e5792..5e9b8f2 100644 --- a/simulate.py +++ b/simulate.py @@ -143,14 +143,18 @@ def costruisci_solutore_implicito_2d( materiale: dict, h_esterno_W_m2K: float, h_interno_W_m2K: float, + h_bordi_W_m2K: float, + lunghezza_conduzione_y_m: float, ): # Costruisce e fattorizza (LU sparsa) la matrice A per Eulero implicito 2D: # A * T_next = rhs # # Le incognite sono i centri cella T[i, j] con i lungo x e j lungo z, - # appiattiti in ordine C (indice = i * n_z + j). I bordi x = 0 e - # x = lunghezza sono adiabatici; i bordi z = 0 (esterno) e z = spessore - # (interno) includono la convezione verso l'ambiente. + # appiattiti in ordine C (indice = i * n_z + j). Tutti e quattro i lati + # della sezione includono la convezione verso l'ambiente; su ogni cella + # agisce inoltre la conduzione circonferenziale (y) verso il resto della + # fascetta, assunto a temperatura ambiente e modellato come + # q_y = -k * (T - T_amb) / L_y². k = materiale["conducibilita_termica_W_mK"] rho = materiale["densita_kg_m3"] cp = materiale["calore_specifico_J_kgK"] @@ -160,21 +164,99 @@ def costruisci_solutore_implicito_2d( r_z = alpha * dt_s / (dz_m * dz_m) b_esterno = h_esterno_W_m2K * dt_s / (rho * cp * dz_m) b_interno = h_interno_W_m2K * dt_s / (rho * cp * dz_m) + b_bordo = h_bordi_W_m2K * dt_s / (rho * cp * dx_m) + c_y = alpha * dt_s / (lunghezza_conduzione_y_m * lunghezza_conduzione_y_m) n = n_x * n_z - convezione = np.zeros(n) - convezione[0::n_z] = b_esterno - convezione[n_z - 1::n_z] = b_interno + scambio = np.full(n, c_y) + scambio[0::n_z] += b_esterno + scambio[n_z - 1::n_z] += b_interno + scambio[:n_z] += b_bordo + scambio[n - n_z:] += b_bordo A = ( sp.identity(n) + r_x * sp.kron(_laplaciano_1d(n_x), sp.identity(n_z)) + r_z * sp.kron(sp.identity(n_x), _laplaciano_1d(n_z)) - + sp.diags(convezione) + + sp.diags(scambio) ) return splu(sp.csc_matrix(A)) +def prepara_stato_termico(fascetta: dict, aria: dict, sorgente: dict) -> dict: + # Prepara griglia, coefficienti e solutore fattorizzato per un run: + # tutto ciò che resta costante durante l'integrazione temporale. + materiale = MATERIALI[fascetta["materiale"]] + + lunghezza = fascetta["lunghezza_mm"] / 1000.0 + spessore = fascetta["spessore_mm"] / 1000.0 + n_x = fascetta["n_nodi_x"] + n_z = fascetta["n_nodi_z"] + dx = lunghezza / n_x + dz = spessore / n_z + lunghezza_y = fascetta["lunghezza_conduzione_y_mm"] / 1000.0 + + dt = SIMULAZIONE["dt_interno_s"] + + 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"]) + + rho = materiale["densita_kg_m3"] + cp = materiale["calore_specifico_J_kgK"] + alpha = materiale["conducibilita_termica_W_mK"] / (rho * cp) + + z_centri = (np.arange(n_z) + 0.5) * dz + + return { + "n_x": n_x, + "n_z": n_z, + "dx_m": dx, + "dz_m": dz, + "dt_s": dt, + "x_centri_m": (np.arange(n_x) + 0.5) * dx, + "z_centri_m": z_centri, + "skin_depth_m": skin_depth, + "rho": rho, + "cp": cp, + "profilo_z": profilo_deposizione_z_1_m(z_centri, spessore, skin_depth), + "b_esterno": aria["h_esterno_W_m2K"] * dt / (rho * cp * dz), + "b_interno": aria["h_interno_W_m2K"] * dt / (rho * cp * dz), + "b_bordo": aria["h_bordi_W_m2K"] * dt / (rho * cp * dx), + "c_y": alpha * dt / (lunghezza_y * lunghezza_y), + "T_ambiente_C": aria["temperatura_ambiente_C"], + "solutore": costruisci_solutore_implicito_2d( + n_x=n_x, + n_z=n_z, + dt_s=dt, + dx_m=dx, + dz_m=dz, + materiale=materiale, + h_esterno_W_m2K=aria["h_esterno_W_m2K"], + h_interno_W_m2K=aria["h_interno_W_m2K"], + h_bordi_W_m2K=aria["h_bordi_W_m2K"], + lunghezza_conduzione_y_m=lunghezza_y, + ), + } + + +def passo_implicito(stato: dict, T: np.ndarray, q_x: np.ndarray) -> np.ndarray: + # Avanza il campo di temperatura di un passo dt: assembla il termine noto + # (sorgente volumetrica, convezione sui quattro lati, conduzione + # circonferenziale verso l'ambiente) e risolve il sistema implicito. + T_amb = stato["T_ambiente_C"] + rhs = T + (stato["dt_s"] / (stato["rho"] * stato["cp"])) * ( + q_x[:, None] * stato["profilo_z"][None, :] + ) + rhs += stato["c_y"] * T_amb + rhs[:, 0] += stato["b_esterno"] * T_amb + rhs[:, -1] += stato["b_interno"] * T_amb + rhs[0, :] += stato["b_bordo"] * T_amb + rhs[-1, :] += stato["b_bordo"] * T_amb + return stato["solutore"].solve(rhs.ravel()).reshape(T.shape) + + def quantizza(valore: float, passo: float) -> float: if passo <= 0.0: return valore @@ -235,49 +317,21 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: sensore = cfg_run["sensore"] nome_materiale = fascetta["materiale"] - materiale = MATERIALI[nome_materiale] - lunghezza = fascetta["lunghezza_mm"] / 1000.0 - spessore = fascetta["spessore_mm"] / 1000.0 - n_x = fascetta["n_nodi_x"] - n_z = fascetta["n_nodi_z"] - dx = lunghezza / n_x - dz = spessore / n_z - x_centri = (np.arange(n_x) + 0.5) * dx - z_centri = (np.arange(n_z) + 0.5) * dz + stato = prepara_stato_termico(fascetta, aria, sorgente) + n_x = stato["n_x"] + n_z = stato["n_z"] + x_centri = stato["x_centri_m"] + skin_depth = stato["skin_depth_m"] x_sensore = sensore["x_mm"] / 1000.0 - i_sensore = min(n_x - 1, max(0, int(x_sensore / dx))) + i_sensore = min(n_x - 1, max(0, int(x_sensore / stato["dx_m"]))) - dt = SIMULAZIONE["dt_interno_s"] + dt = stato["dt_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"]) - - solutore = costruisci_solutore_implicito_2d( - n_x=n_x, - n_z=n_z, - dt_s=dt, - dx_m=dx, - dz_m=dz, - materiale=materiale, - h_esterno_W_m2K=aria["h_esterno_W_m2K"], - h_interno_W_m2K=aria["h_interno_W_m2K"], - ) - - rho = materiale["densita_kg_m3"] - cp = materiale["calore_specifico_J_kgK"] - - b_esterno = aria["h_esterno_W_m2K"] * dt / (rho * cp * dz) - b_interno = aria["h_interno_W_m2K"] * dt / (rho * cp * dz) - - profilo_z = profilo_deposizione_z_1_m(z_centri, spessore, skin_depth) - - T = np.full((n_x, n_z), fascetta["temperatura_iniziale_C"], dtype=float) + T = np.full((n_x, n_z), aria["temperatura_ambiente_C"], dtype=float) T_sensore = T[i_sensore, -1] prossimo_campione_t = 0.0 @@ -310,12 +364,7 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: x_sorgente, q_x = profilo_flusso_incidente_W_m2( sorgente, x_sensore, t, x_centri ) - - rhs = T + (dt / (rho * cp)) * q_x[:, None] * profilo_z[None, :] - rhs[:, 0] += b_esterno * aria["temperatura_ambiente_C"] - rhs[:, -1] += b_interno * aria["temperatura_ambiente_C"] - - T = solutore.solve(rhs.ravel()).reshape(n_x, n_z) + T = passo_implicito(stato, T, q_x) # Temperatura vera della superficie interna nel punto osservato # dal sensore infrarosso. @@ -358,8 +407,9 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: "file_csv": str(output_csv.name), "materiale": nome_materiale, "diametro_m": fascetta["diametro_mm"] / 1000.0, - "lunghezza_m": lunghezza, - "spessore_m": spessore, + "lunghezza_m": fascetta["lunghezza_mm"] / 1000.0, + "spessore_m": fascetta["spessore_mm"] / 1000.0, + "lunghezza_conduzione_y_m": fascetta["lunghezza_conduzione_y_mm"] / 1000.0, "n_nodi_x": n_x, "n_nodi_z": n_z, "durata_s": durata, @@ -368,6 +418,7 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: "temperatura_ambiente_C": aria["temperatura_ambiente_C"], "h_esterno_W_m2K": aria["h_esterno_W_m2K"], "h_interno_W_m2K": aria["h_interno_W_m2K"], + "h_bordi_W_m2K": aria["h_bordi_W_m2K"], "x_inizio_m": sorgente["x_inizio_m"], "x_fine_m": sorgente["x_fine_m"], "x_sensore_m": x_sensore,