diff --git a/.gitignore b/.gitignore index fc5ba3c..e5d24d3 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,3 @@ .venv/ __pycache__/ -dataset/ +output/ diff --git a/CLAUDE.md b/CLAUDE.md index ed2a9e0..31c8ba4 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -20,14 +20,10 @@ pip install -r requirements.txt # Punto di ingresso unico: `python main.py` elenca le azioni disponibili python main.py mesh # disegna la sola mesh a elementi shell -python main.py simula # genera dataset/run_XXXX.csv + dataset/metadata.csv -python main.py grafico # grafici del primo run -python main.py anima # animazione 2D della sezione -python main.py anima3d # animazione 3D isometrica del barattolo python main.py fem # animazione 3D del campo FEM sulla mesh shell ``` -Ogni modulo resta eseguibile anche direttamente (`python simulate.py`, `python plot_csv.py`, ...). +Ogni modulo resta eseguibile anche direttamente (`python plot_mesh.py`, `python plot_animazione_fem.py`). Attivare sempre il venv (`source .venv/bin/activate`) prima di eseguire qualsiasi comando Python. @@ -35,52 +31,23 @@ Non sono configurati test o linter. ## Architettura -Generatore di dataset per misurazioni termiche pseudo-realistiche di una fascetta (anello cilindrico sottile) riscaldata da sorgenti a induzione in movimento. +Analisi termica transitoria a elementi finiti di una fascetta (anello cilindrico sottile) riscaldata da sorgenti a induzione in movimento. -**Geometria:** la fascetta ha diametro, spessore e lunghezza configurabili. Il dominio simulato è la sezione rettangolare lunghezza × spessore, con origine (0, 0) nel vertice in alto a sinistra: x = lunghezza (le sorgenti si muovono in direzione -x sul lato esterno), z = spessore (0 = lato esterno, spessore = lato interno). La coordinata circonferenziale y non è risolta: l'offset y delle sorgenti è collassato in un'attenuazione gaussiana del flusso. Il sensore è un pirometro a infrarossi dentro la fascetta che misura la superficie interna in un punto x fisso. +**Geometria:** la fascetta ha diametro, spessore e lunghezza configurabili. Il dominio discretizzato è la superficie media del cilindro: x = asse (le sorgenti si muovono lungo x sulla superficie esterna), theta = coordinata circonferenziale con theta = 0 sul piano del sensore, s = R·theta la lunghezza d'arco. Lo spessore non è discretizzato. Il sensore è un pirometro a infrarossi dentro la fascetta che misura la parete in un punto x fisso sul piano theta = 0. **Flusso dei dati:** -0. `main.py` — dispatcher da riga di comando che seleziona l'azione (mesh, simula, grafico, anima, anima3d) -1. `config.py` — tutti i parametri configurabili (dizionari SIMULAZIONE, FASCETTA, MESH, ARIA, SORGENTE, SENSORE, RANDOMIZZAZIONE) -2. `materials.py` — dizionario MATERIALI con proprietà termofisiche ed elettriche per materiale -3. `simulate.py` — motore principale: genera N run randomizzati, scrive i CSV, scrive `metadata.csv` -4. `plot_csv.py` — visualizzazione autonoma per un singolo run -5. `mesh.py` — mesh a elementi shell quadrangolari della superficie media del cilindro, visualizzata da `plot_mesh.py` -6. `fem.py` — solutore termico transitorio a elementi finiti sulla mesh shell, animato da `plot_animazione_fem.py` +1. `main.py` — dispatcher da riga di comando che seleziona l'azione (mesh, fem) +2. `config.py` — tutti i parametri configurabili (dizionari FASCETTA, MESH, FEM, ARIA, SORGENTE, SENSORE, RANDOMIZZAZIONE, USCITA) +3. `materials.py` — dizionario MATERIALI con proprietà termofisiche ed elettriche per materiale +4. `mesh.py` — mesh a elementi shell quadrangolari della superficie media del cilindro, visualizzata da `plot_mesh.py` +5. `fem.py` — solutore termico transitorio a elementi finiti sulla mesh shell, animato da `plot_animazione_fem.py` -**Mesh shell (`mesh.py`):** griglia strutturata `n_elementi_x × n_elementi_circonferenza` di quadrilateri a 4 nodi sulla superficie media (raggio = (diametro − spessore)/2); lo spessore non è discretizzato, è un attributo della shell. La mesh è chiusa lungo la circonferenza (nessun nodo duplicato sulla cucitura) e usa lo stesso sistema di coordinate globale di `plot_animazione_3d.py`: x = asse, y = R·sin(theta), z = R·cos(theta), con theta = 0 sul piano del sensore. +**Mesh shell (`mesh.py`):** griglia strutturata `n_elementi_x × n_elementi_circonferenza` di quadrilateri a 4 nodi sulla superficie media (raggio = (diametro − spessore)/2); lo spessore non è discretizzato, è un attributo della shell. La mesh è chiusa lungo la circonferenza (nessun nodo duplicato sulla cucitura) e usa coordinate globali x = asse, y = R·sin(theta), z = R·cos(theta), con theta = 0 sul piano del sensore. -## Due modelli termici distinti +**Pipeline FEM dentro `fem.py`:** 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 `h_bordi` sui due bordi anulari (area = spessore × perimetro). Il collasso in spessore è lecito perché Biot ≈ 1e-7 e il tempo di diffusione nella parete è di pochi ms: la profondità di penetrazione del riscaldamento non entra nel modello, conta solo il flusso assorbito. La skin depth è comunque calcolata dalla frequenza dell'induttore (`calcola_skin_depth_m`) e riportata in `stato["skin_depth_m"]` come sola diagnostica. L'impronta della sorgente è una gaussiana isotropa in (x, s) con distanza circonferenziale valutata sull'immagine più vicina (superficie chiusa). Gli elementi sono bilineari a 4 nodi e tutti identici (rettangoli dx × ds), quindi le matrici di elemento (`MASSA_RIF`, `RIGIDEZZA_X`, `RIGIDEZZA_S`) sono calcolate una volta e assemblate in forma vettorizzata; la matrice di massa di superficie è riusata per capacità, convezione e carico della sorgente. Eulero implicito con matrice costante fattorizzata LU una volta per analisi (`prepara_stato_fem`), poi `passo_implicito_fem` assembla solo il termine noto. In `fem.py` stanno anche la cinematica del gruppo di sorgenti (`_x_riferimento_iniziale_m`, `_x_riferimento_finale_m`, `_intervallo_attivo`) e la lettura della configurazione (`configurazione_randomizzata`). -Il progetto contiene due discretizzazioni indipendenti della stessa fisica; scambiano solo la cinematica delle sorgenti (`_x_riferimento_iniziale_m`, `_x_riferimento_finale_m`, `_intervallo_attivo` in [simulate.py](simulate.py)), che deve restare unica. - -| | volumi finiti ([simulate.py](simulate.py)) | elementi finiti ([fem.py](fem.py)) | -|---|---|---| -| dominio | sezione x-z (spessore risolto) | superficie media x-s (circonferenza risolta) | -| circonferenza y | attenuazione gaussiana + sink di aletta | direzione risolta, conduzione reale | -| spessore | `n_nodi_z` celle, deposizione esponenziale da skin depth | collassato: T uniforme nella parete | -| griglia | `FASCETTA["n_nodi_x"]`, `n_nodi_z` | `MESH["n_elementi_x"]`, `n_elementi_circonferenza` | -| passo temporale | `SIMULAZIONE["dt_interno_s"]` | `FEM["dt_s"]` | -| output | dataset CSV + animazioni 2D/3D | animazione 3D sulla mesh | - -**Pipeline FEM dentro `fem.py`:** 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 `h_bordi` sui due bordi anulari (area = spessore × perimetro). Il collasso in spessore è lecito perché Biot ≈ 1e-7 e il tempo di diffusione nella parete è di pochi ms: la skin depth non entra nel modello shell, conta solo il flusso assorbito. L'impronta della sorgente è una gaussiana isotropa in (x, s) con distanza circonferenziale valutata sull'immagine più vicina (superficie chiusa). Gli elementi sono bilineari a 4 nodi e tutti identici (rettangoli dx × ds), quindi le matrici di elemento (`MASSA_RIF`, `RIGIDEZZA_X`, `RIGIDEZZA_S`) sono calcolate una volta e assemblate in forma vettorizzata; la matrice di massa di superficie è riusata per capacità, convezione e carico della sorgente. Eulero implicito con matrice costante fattorizzata LU una volta per run (`prepara_stato_fem`), poi `passo_implicito_fem` assembla solo il termine noto. - -**Pipeline fisica dentro `simula_singolo()` in [simulate.py](simulate.py):** - -- 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: `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. Il termine è un'equazione di aletta ricavata sull'intero volume del cilindro: il calore conduce lungo y attraverso l'intero spessore mentre le superfici esterna e interna dell'intero cilindro scambiano per convezione, dando q_y = -(h_esterno + h_interno)/spessore · (T - T_amb), senza parametri di conduzione y configurabili a parte -- 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. +**Sensore e randomizzazione:** il campo nodale è la temperatura vera della parete; `simula_campo_fem` restituisce anche `T_lette`, la lettura del sensore degradata da inerzia del primo ordine (a ogni passo dt), rumore gaussiano e quantizzazione (a ogni campionamento di frame). `configurazione_randomizzata` perturba velocità, flusso di picco, sigma dello spot, offset y, temperatura ambiente e rumore del sensore con un RNG con seed `FEM["seed"]`, così due analisi restano riproducibili ma diverse; con `RANDOMIZZAZIONE["abilitata"] = False` restituisce i valori nominali di `config.py`. ## Convenzioni su `config.py` @@ -88,8 +55,8 @@ Ogni parametro in [config.py](config.py) ha un commento che spiega solo cos'è ( ## Vincoli progettuali chiave -- Il modello ai volumi finiti di `simulate.py` è 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. Il modello FEM di `fem.py` fa il contrario (circonferenza risolta, spessore collassato): non modificare uno assumendo che valgano le ipotesi dell'altro. +- Il modello è una shell: la circonferenza è risolta, lo spessore è collassato (temperatura uniforme nella parete). Non introdurre effetti attraverso lo spessore (skin depth, gradiente radiale) senza prima discretizzarlo. - 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`). +- La matrice implicita è costruita e fattorizzata una volta sola (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. +- Le visualizzazioni salvano in `USCITA["cartella"]` solo quando il backend matplotlib non è interattivo; altrimenti aprono una finestra. - Aggiungere un nuovo materiale richiede solo una nuova voce nel dizionario `MATERIALI` in [materials.py](materials.py); la chiave del materiale va poi impostata in `FASCETTA["materiale"]` in [config.py](config.py). diff --git a/README.md b/README.md index a0c7e71..79f6940 100644 --- a/README.md +++ b/README.md @@ -1,113 +1,132 @@ -# Simulatore Termico 2D — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso +# Analisi Termica FEM — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso -Questo progetto genera misurazioni CSV pseudo-realistiche della temperatura di una -fascetta (anello cilindrico sottile) riscaldata da un gruppo di sorgenti a induzione -in movimento, osservata da un sensore a infrarossi fisso. Lo scopo è produrre dataset -per l'addestramento e la validazione di modelli di stima/regressione termica. +Questo progetto calcola il transitorio termico di una fascetta (anello cilindrico +sottile) riscaldata da un gruppo di sorgenti a induzione in movimento, osservata da un +sensore a infrarossi fisso. La fascetta è discretizzata con una mesh a elementi shell +quadrangolari e il campo di temperatura è risolto a elementi finiti, poi visualizzato in +un'animazione 3D isometrica. ## Geometria La fascetta è un anello cilindrico definito da tre dimensioni: -- **diametro** (default 70 mm) — il diametro del cilindro; -- **spessore** (default 0.12 mm) — lo spessore della parete; +- **diametro** (default 70 mm) — il diametro esterno del cilindro; +- **spessore** (default 0.18 mm) — lo spessore della parete; - **lunghezza** (default 100 mm) — l'estensione assiale lungo `x`. -Il dominio simulato è la **sezione rettangolare lunghezza × spessore**. Il sistema di -coordinate ha l'origine `(0, 0)` nel vertice in alto a sinistra della sezione: +Il dominio discretizzato è la **superficie media** del cilindro (raggio = +(diametro − spessore)/2). Lo spessore non è discretizzato: è un attributo degli +elementi shell. ```text - sorgenti (induttori), in moto verso -x + sorgenti (induttori), in moto lungo x sulla superficie esterna ▼ ▼ ▼ - (0,0) ─────────────────────────────────────► x - │ ┌───────────────────────────────────┐ z = 0 lato ESTERNO - │ │ sezione della fascetta │ (flusso termico) - │ └───────────────────────────────────┘ z = spessore lato INTERNO - ▼ ┆ - z ┆ linea di vista - ▲ - sensore IR (fisso, x = 50 mm, - a 10 mm dalla parete interna) + ┌───────────────────────────────────┐ + x = 0 │ superficie media del cilindro │ x = lunghezza + └───────────────────────────────────┘ + ┆ linea di vista + ▲ + sensore IR (fisso, x = 50 mm, theta = 0, + a 10 mm dalla parete interna) ``` -- **x** = direzione della lunghezza, da `0` a `lunghezza`. Le sorgenti viaggiano in - direzione `-x` sul lato esterno. -- **z** = direzione dello spessore, da `0` (lato esterno, dove arriva il flusso - termico) a `spessore` (lato interno, osservato dal sensore). -- **y** = coordinata circonferenziale (lungo la circonferenza π·diametro). Non è - risolta spazialmente: vedi sotto come viene trattata. +- **x** = asse del cilindro, da `0` a `lunghezza`. +- **theta** = coordinata circonferenziale, con `theta = 0` sul piano del sensore; + nelle formule si usa la lunghezza d'arco `s = R·theta`. +- Coordinate globali per la visualizzazione: `x`, `y = R·sin(theta)`, + `z = R·cos(theta)`. Il **sensore** è un pirometro a infrarossi posto all'interno della fascetta, a una -distanza configurabile dalla parete interna (default 10 mm). Essendo senza contatto, -la distanza non influenza la misura: il sensore legge la temperatura della superficie -interna nel punto `x` configurato (default 50 mm, al centro della lunghezza). +distanza configurabile dalla parete interna (default 10 mm). Essendo senza contatto, la +distanza non influenza la misura: il sensore legge la temperatura della parete nel nodo +più vicino a `x = 50 mm` sul piano `theta = 0`. ## Modello fisico -Non è una simulazione FEM elettromagnetica + termica completa: è un generatore -pratico di dataset. La catena di approssimazioni è la seguente. +L'equazione risolta è quella del calore per unità di superficie media: -1. **Sorgenti gaussiane in moto** — ogni sorgente ha un'impronta gaussiana di raggio - `sigma_punto_m`. Un gruppo di `numero_sorgenti` sorgenti equidistanti - (`distanza_sorgenti_m`) si muove rigidamente a velocità costante. Il profilo di - flusso sul lato esterno è la somma dei contributi: - `q(x, t) = Σᵢ q_picco · efficienza · exp(-((x - xᵢ(t))² + Δy²) / (2σ²))`. - L'offset circonferenziale `Δy` tra il percorso delle sorgenti e il punto osservato - dal sensore non è risolto spazialmente: entra come attenuazione gaussiana del flusso. +```text +rho·cp·t ∂T/∂t = k·t ∇²T + q(x, s, t) − (h_est + h_int)·(T − T_amb) +``` -2. **Skin depth** — il riscaldamento a induzione è approssimato come riscaldamento - volumetrico che decade esponenzialmente con la profondità `z`: - `q_vol(x, z) = q(x) · exp(-z/δ) / (δ·(1 - exp(-spessore/δ)))`, normalizzato in modo - da conservare il flusso superficiale. La skin depth `δ = √(2ρₑ/(ωμ))` è calcolata - dalla resistività elettrica e dalla permeabilità del materiale alla frequenza di - induzione, oppure può essere imposta con `skin_depth_fissa_m`. Per la banda - stagnata a 20 kHz risulta ≈ 0.1 mm, confrontabile con lo spessore: la parete è - quasi isoterma attraverso lo spessore. +con `t` = spessore della parete e `∇²` il laplaciano nel piano `(x, s)`. -3. **Diffusione 2D del calore** — l'equazione del calore è integrata nella sezione - `(x, z)` con volumi finiti ed Eulero implicito (incondizionatamente stabile). +1. **Collasso dello spessore** — la temperatura è uniforme attraverso la parete. 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 + `δ = √(2ρₑ/(ωμ))` è comunque calcolata dalla frequenza dell'induttore + (`frequenza_hz`, oppure imposta con `skin_depth_fissa_m`) e riportata come + grandezza diagnostica: per la banda stagnata a 20 kHz risulta ≈ 0.14 mm, + confrontabile con lo spessore, il che conferma che la parete è quasi isoterma + attraverso lo spessore. -4. **Scambi con l'esterno** — la sezione scambia calore con l'ambiente su tutto il - contorno: - - convezione sul lato esterno (`h_esterno_W_m2K`), sul lato interno - (`h_interno_W_m2K`) e sui due bordi in x (`h_bordi_W_m2K`); - - **conduzione circonferenziale**: la sezione cede calore per conduzione lungo `y` - al resto della fascetta, assunto a temperatura ambiente. Il termine è - un'equazione di aletta ricavata sull'intero volume del cilindro: il calore - conduce lungo `y` attraverso l'intero spessore mentre le superfici esterna e - interna dell'intero cilindro perdono calore per convezione, dando - `q_y = -(h_esterno + h_interno)/spessore · (T - T_amb)` — nessun parametro di - conduzione `y` aggiuntivo da configurare. +2. **Sorgenti gaussiane in moto** — ogni sorgente ha un'impronta gaussiana isotropa nel + piano `(x, s)` di raggio `sigma_punto_m`. Un gruppo di `numero_sorgenti` sorgenti + equidistanti (`distanza_sorgenti_m`) si muove rigidamente a velocità costante: + `q(x, s, t) = Σᵢ q_picco · efficienza · exp(-((x - xᵢ(t))² + (s - s₀)²) / (2σ²))`. + La distanza circonferenziale è valutata sull'immagine più vicina, perché la + superficie è chiusa. L'offset `s₀` del percorso rispetto al piano del sensore è + `offset_y_percorso_m`. -5. **Temperatura iniziale** — il campo parte uniformemente alla temperatura ambiente - del run (che è randomizzata, quindi varia run per run). +3. **Scambi con l'ambiente** — convezione sulla faccia esterna (`h_esterno_W_m2K`) e su + quella interna (`h_interno_W_m2K`) su tutta la superficie, più convezione sui due + bordi anulari `x = 0` e `x = lunghezza` (`h_bordi_W_m2K`, su un'area pari a + spessore × perimetro). -6. **Sensore realistico** — la lettura aggiunge alla temperatura vera della superficie - interna: inerzia del primo ordine (`costante_tempo_s`), rumore gaussiano - (`rumore_std_C`) e quantizzazione (`quantizzazione_C`). +4. **Temperatura iniziale** — il campo parte uniformemente alla temperatura ambiente + dell'analisi (randomizzabile). + +5. **Sensore reale** — la lettura aggiunge alla temperatura vera della parete nel punto + osservato: inerzia del primo ordine (`costante_tempo_s`), rumore gaussiano + (`rumore_std_C`) e quantizzazione (`quantizzazione_C`), gli ultimi due applicati a + ogni campionamento. + +6. **Randomizzazione** — con `RANDOMIZZAZIONE["abilitata"]` ogni analisi perturba + velocità, flusso di picco, sigma dello spot, offset y del percorso, temperatura + ambiente e rumore del sensore con estrazioni gaussiane/uniformi da un RNG con seed + fisso (`FEM["seed"]`): analisi diverse dello stesso setup differiscono come + differirebbero due passaggi reali, restando riproducibili. + +Irraggiamento non modellato e proprietà dei materiali costanti con la temperatura: a +~210 °C le perdite radiative non sono del tutto trascurabili rispetto alla convezione, e +per gli acciai ferromagnetici la skin depth reale varia fortemente con temperatura e +campo (punto di Curie non modellato). ## Metodo numerico -- Griglia a volumi finiti `n_nodi_x × n_nodi_z` (default 100 × 15); le incognite sono - i centri cella. -- Eulero implicito con passo `dt_interno_s` (default 0.2 ms), più fine del periodo di - campionamento CSV. -- Tutti i termini (diffusione, convezione, conduzione circonferenziale) sono lineari e - costanti nel run: la matrice sparsa viene costruita e **fattorizzata LU una sola - volta per run** (`scipy.sparse.linalg.splu`); ogni passo temporale risolve solo i - sistemi triangolari. Un run da 30 s simulati richiede ~20 s di calcolo. +- **Mesh** (`mesh.py`): griglia strutturata `n_elementi_x × n_elementi_circonferenza` + di quadrilateri a 4 nodi (default 40 × 48), chiusa lungo la circonferenza senza nodi + duplicati sulla cucitura. +- **Elementi**: shell bilineari a 4 nodi. Tutti gli elementi sono rettangoli identici + `dx × ds`, quindi le matrici di elemento (massa e rigidezza) sono calcolate una volta + sola e assemblate in forma vettorizzata. La matrice di massa di superficie è riusata + per capacità termica, convezione sulle facce e carico della sorgente; la convezione + sui bordi anulari usa una matrice di massa 1D sull'anello di nodi. +- **Integrazione temporale**: Eulero implicito con passo `FEM["dt_s"]` (default 1 ms). + Tutti i termini sono lineari e costanti, quindi la matrice di sistema è **assemblata e + fattorizzata LU una sola volta** (`scipy.sparse.linalg.splu`) in `prepara_stato_fem`; + ogni passo assembla solo il termine noto e risolve i sistemi triangolari. +- **Flusso nodale**: l'impronta gaussiana isotropa è separabile, + `q(x, s) = q_x(x) · q_s(s)`, e il fattore circonferenziale è costante nel tempo + (precalcolato). A ogni passo si valuta quindi solo il fattore assiale sui `n_x + 1` + nodi distinti in x, espanso con un prodotto esterno. - Se in futuro si introducessero proprietà dipendenti dalla temperatura, la matrice andrebbe ricostruita e rifattorizzata a ogni passo. ## File ```text -config.py tutti i parametri di simulazione -materials.py proprietà termofisiche ed elettriche dei materiali -simulate.py motore fisico + generazione dei CSV -plot_csv.py grafici rapidi (temperature e flusso) del primo run -plot_animazione.py animazione della sezione: campo T(x,z), sorgenti, sensore -dataset/ output generato da simulate.py (ricreato a ogni esecuzione) +main.py punto di ingresso da riga di comando +config.py tutti i parametri (geometria, mesh, FEM, aria, sorgente, sensore) +materials.py proprietà termofisiche ed elettriche dei materiali +mesh.py generazione della mesh a elementi shell +fem.py solutore termico transitorio a elementi finiti +plot_mesh.py disegno della sola mesh +plot_animazione_fem.py animazione 3D del campo di temperatura +output/ immagini e GIF salvate quando il backend non è interattivo ``` ## Installazione @@ -121,38 +140,38 @@ pip install -r requirements.txt ## Uso ```bash -# genera il dataset (ATTENZIONE: cancella e ricrea la cartella dataset/) -python simulate.py - -# grafici statici del primo run (temperatura e flusso nel tempo) -python plot_csv.py - -# animazione della sezione durante il passaggio delle sorgenti -python plot_animazione.py +python main.py # elenco delle azioni disponibili +python main.py mesh # disegna la sola mesh a elementi shell +python main.py fem # integra il campo FEM e lo anima in 3D ``` -Gli script di visualizzazione aprono finestre interattive (backend Qt); se il backend -non è interattivo (es. sessione senza display) salvano automaticamente PNG/GIF in -`dataset/`. +Ogni modulo resta eseguibile anche direttamente (`python plot_mesh.py`, +`python plot_animazione_fem.py`). -L'animazione riproduce esattamente la fisica di `run_0001` (stesso seed) e mostra tre -pannelli allineati: il profilo di flusso `q(x)` con le sorgenti in transito, il campo -di temperatura nella sezione con il sensore IR, e la temperatura nel punto osservato -(vera e con inerzia del sensore). Finestra temporale e cadenza dei fotogrammi si -regolano con le costanti in testa a `plot_animazione.py`. +Gli script aprono finestre interattive (backend Qt); se il backend non è interattivo +(es. sessione senza display) salvano automaticamente PNG/GIF in `output/`. + +L'animazione FEM mostra due pannelli: la fascetta in vista isometrica con il campo di +temperatura dipinto sugli elementi shell e i marker delle sorgenti in transito, e la +storia della temperatura nel punto osservato dal sensore (valore nodale vero e lettura +del sensore reale). Finestra temporale iniziale e cadenza dei fotogrammi si +regolano con le costanti in testa a `plot_animazione_fem.py`; la durata simulata è +`FEM["durata_s"]`. ## Configurazione Tutto si modifica in `config.py`. I dizionari principali: -| Dizionario | Contenuto | -|------------------|---------------------------------------------------------------------------| -| `SIMULAZIONE` | numero di run, durata, campionamento CSV, passo interno, seed, cartella | -| `FASCETTA` | diametro, lunghezza, spessore, griglia, conduzione circonferenziale, materiale | -| `ARIA` | temperatura ambiente e coefficienti di convezione dei quattro lati | -| `SORGENTE` | corsa, velocità (il segno dà il verso), gruppo di sorgenti, gaussiana, flusso, frequenza | -| `SENSORE` | posizione (x e distanza dalla parete), inerzia, rumore, quantizzazione | -| `RANDOMIZZAZIONE`| entità delle perturbazioni per run | +| Dizionario | Contenuto | +|-------------|---------------------------------------------------------------------------| +| `FASCETTA` | diametro, lunghezza, spessore, materiale | +| `MESH` | numero di elementi shell lungo x e lungo la circonferenza | +| `FEM` | passo temporale, durata dell'analisi transitoria, seed | +| `ARIA` | temperatura ambiente e coefficienti di convezione (facce e bordi) | +| `SORGENTE` | corsa, velocità (il segno dà il verso), gruppo di sorgenti, gaussiana, flusso, frequenza | +| `SENSORE` | posizione (x e distanza dalla parete), inerzia, rumore, quantizzazione | +| `RANDOMIZZAZIONE` | entità delle perturbazioni per analisi | +| `USCITA` | cartella per immagini e animazioni salvate | Punti da conoscere: @@ -163,54 +182,14 @@ Punti da conoscere: alla fine della propria corsa. - **Materiale**: `FASCETTA["materiale"]` deve essere una chiave di `MATERIALI` in `materials.py`. Per aggiungere un materiale basta una nuova voce nel dizionario - (conducibilità, densità, calore specifico, resistività elettrica, permeabilità). + (conducibilità termica, densità, calore specifico, resistività elettrica, + permeabilità relativa). - **Proprietà di `banda_stagnata`**: la banda stagnata è un nastro di acciaio a basso tenore di carbonio (0,15–0,25% C) rivestito su entrambe le facce da un sottile strato di stagno elettrolitico, dello spessore di pochi micrometri — trascurabile rispetto allo spessore tipico della parete (es. 0,18 mm) e quindi ininfluente sulle - proprietà termiche, elettriche e magnetiche in massa. I valori in `materials.py` - sono quindi quelli dell'acciaio dolce sottostante: + proprietà termiche, elettriche e magnetiche in massa. I valori in `materials.py` sono + quindi quelli dell'acciaio dolce sottostante: [Banda stagnata: tutto quello che c'è da sapere (MUNDOLATAS)](https://mundolatas.com/it/banda-stagnata-tutto-quello-che-ce-da-sapere/), [Bande stagnate elettrolitiche (EUROPERF)](https://www.europerf.it/it/banda-stagnata-elettrolitica.php), [Differenza tra banda stagnata e acciaio inossidabile (Wuxi Bright Packing)](https://it.brightmetalcan.com/info/difference-between-tinplate-and-stainless-stee-48700260.html). -- **Randomizzazione**: ogni run perturba velocità, flusso di picco, sigma, offset y, - temperatura ambiente e rumore del sensore con estrazioni da un RNG a seed fisso - (`SIMULAZIONE["seed"]`): il dataset è riproducibile. - -## Output - -### `dataset/run_XXXX.csv` — serie temporale del run - -| Colonna | Significato | -|--------------------------------|--------------------------------------------------------------------| -| `id_run` | identificativo del run | -| `tempo_s` | tempo simulato | -| `x_sorgente_m` | posizione della sorgente di riferimento del gruppo | -| `offset_y_sorgente_m` | offset circonferenziale del percorso (costante nel run) | -| `flusso_termico_sorgente_W_m2` | flusso efficace nel punto x del sensore | -| `skin_depth_m` | skin depth usata (costante nel run) | -| `T_vera_lato_sensore_C` | temperatura vera della superficie interna nel punto del sensore | -| `T_misurata_sensore_C` | lettura del sensore (inerzia + rumore + quantizzazione) | -| `T_lato_caldo_C` | temperatura della superficie esterna nello stesso punto x | -| `T_ambiente_C` | temperatura ambiente del run | -| `velocita_m_s`, `sigma_punto_m`, `flusso_picco_W_m2` | parametri randomizzati del run | -| `materiale` | chiave del materiale | - -### `dataset/metadata.csv` — una riga per run - -Contiene tutti i parametri effettivi del run (geometria, griglia, coefficienti di -scambio, parametri delle sorgenti e del sensore, valori randomizzati) e le temperature -di picco vera e misurata: utile come ground truth e per filtrare i run. - -## Limitazioni - -1. Il campo elettromagnetico non è simulato: l'accoppiamento induttivo è ridotto a - impronta gaussiana × efficienza × decadimento esponenziale in z. -2. La coordinata circonferenziale y non è risolta: offset del percorso e conduzione - verso il resto della fascetta sono modelli collassati (attenuazione gaussiana e - scambio lineare verso T ambiente). -3. Le proprietà dei materiali sono costanti con la temperatura; per gli acciai - ferromagnetici la skin depth reale varia fortemente con temperatura e campo - (punto di Curie non modellato). -4. Irraggiamento non modellato: a ~220 °C le perdite radiative non sono del tutto - trascurabili rispetto alla convezione. diff --git a/config.py b/config.py index f3b6b55..01ccaba 100644 --- a/config.py +++ b/config.py @@ -1,30 +1,21 @@ -# Configurazione per il simulatore termico 2D della sezione di una fascetta. +# Configurazione dell'analisi termica a elementi finiti della fascetta. # # Geometria e modello fisico: # - La fascetta è un anello cilindrico con diametro "diametro_mm", spessore # "spessore_mm" e lunghezza "lunghezza_mm". -# - Il dominio simulato è la sezione rettangolare lunghezza × spessore. -# - Sistema di coordinate: origine (0, 0) nel vertice in alto a sinistra -# della sezione. x = direzione della lunghezza (da 0 a lunghezza), -# z = direzione dello spessore (0 = lato esterno, dove agiscono le -# sorgenti; spessore = lato interno, osservato dal sensore). -# - y è la coordinata circonferenziale: non è risolta spazialmente, l'offset -# y del percorso delle sorgenti è collassato in un'attenuazione gaussiana -# del flusso. -# - 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). Scambia inoltre per conduzione -# lungo y con il resto della fascetta, assunto a temperatura ambiente: -# il calore conduce attraverso l'intero volume dello spessore mentre le -# superfici esterna e interna dell'intero cilindro perdono calore per -# convezione (equazione dell'aletta), derivato da h_esterno, h_interno e -# spessore_mm senza parametri di conduzione y aggiuntivi. +# - Il dominio discretizzato è la superficie media del cilindro (mesh a +# elementi shell quadrangolari, vedi mesh.py): x = asse della fascetta, +# s = R·theta = coordinata circonferenziale. Lo spessore non è +# discretizzato, la temperatura è uniforme attraverso la parete. +# - Le sorgenti a induzione si muovono lungo x sulla superficie esterna, con +# un'impronta gaussiana isotropa nel piano (x, s). +# - La superficie scambia per convezione con l'aria sulla faccia esterna e su +# quella interna, più i due bordi anulari x = 0 e x = lunghezza. # - 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". +# contatto la temperatura della superficie interna nel punto x = "x_mm", +# sul piano circonferenziale theta = 0. # # Unità di misura: # - lunghezza: m (mm dove indicato dal suffisso) @@ -33,32 +24,6 @@ # - flusso termico: W/m² # - coefficiente di convezione: W/(m² K) -SIMULAZIONE = { - # Numero di file CSV da generare. - "num_run": 1, - - # Tempo simulato totale. - "durata_s": 30.0, - - # Frequenza di campionamento CSV. - # Esempio: 2 Hz significa una riga ogni 0.5 s. - "frequenza_campionamento_hz": 10.0, - - # Passo di integrazione numerica interna. - # Può essere inferiore al periodo di campionamento CSV. - "dt_interno_s": 0.0002, - - # Seed per la riproducibilità. - "seed": 42, - - # Cartella di output. - "cartella_output": "dataset", - - # Numero di processi paralleli per la generazione dei run. - # None = usa tutti i core disponibili. - "num_processi": None, -} - FASCETTA = { # Diametro della fascetta [mm]. "diametro_mm": 70.0, @@ -69,13 +34,6 @@ FASCETTA = { # Spessore della parete [mm]. "spessore_mm": 0.18, - # Numero di celle del volume finito lungo x (lunghezza). - "n_nodi_x": 100, - - # Numero di celle del volume finito lungo z (spessore). - # Più nodi = maggiore risoluzione spaziale, simulazione più lenta. - "n_nodi_z": 15, - # Deve corrispondere a una chiave in materials.py. "materiale": "banda_stagnata", } @@ -92,19 +50,25 @@ MESH = { FEM = { # Passo di integrazione temporale del solutore a elementi finiti. "dt_s": 0.001, + + # Tempo simulato totale dell'analisi transitoria. + "durata_s": 30.0, + + # Seed per la riproducibilità di randomizzazione e rumore del sensore. + "seed": 42, } ARIA = { # Temperatura dell'aria ambiente. "temperatura_ambiente_C": 25.0, - # Coefficiente di convezione sul lato esterno (z = 0, lato sorgenti). + # Coefficiente di convezione sulla faccia esterna (lato sorgenti). "h_esterno_W_m2K": 12.0, - # Coefficiente di convezione sul lato interno (z = spessore, lato sensore). + # Coefficiente di convezione sulla faccia interna (lato sensore). "h_interno_W_m2K": 8.0, - # Coefficiente di convezione sui bordi laterali (x = 0 e x = lunghezza). + # Coefficiente di convezione sui bordi anulari (x = 0 e x = lunghezza). "h_bordi_W_m2K": 10.0, } @@ -141,7 +105,8 @@ SORGENTE = { # Frazione del flusso incidente che diventa effettivamente calore nella fascetta. "efficienza_riscaldamento": 0.35, - # Frequenza di induzione usata per stimare la skin depth se skin_depth_fissa_m è None. + # Frequenza di induzione dell'induttore, usata per stimare la skin depth + # se skin_depth_fissa_m è None. "frequenza_hz": 20000.0, # Override della skin depth. Usare None per calcolarla dalle proprietà elettriche del materiale. @@ -156,7 +121,7 @@ SENSORE = { # Coordinata x del punto della superficie interna osservato dal sensore [mm]. "x_mm": 50.0, - # Distanza del sensore dalla parete interna lungo z [mm]. + # Distanza del sensore dalla parete interna lungo lo spessore [mm]. # Il sensore è a infrarossi: la distanza non influenza la misura, # è registrata solo come geometria del setup. "distanza_parete_mm": 10.0, @@ -175,7 +140,7 @@ SENSORE = { } RANDOMIZZAZIONE = { - # Se abilitata, ogni run varia leggermente alcuni parametri. + # Se abilitata, ogni analisi varia leggermente alcuni parametri. "abilitata": False, # Deviazioni standard relative. @@ -192,3 +157,9 @@ RANDOMIZZAZIONE = { # dalla linea ideale allineata con il sensore. "offset_y_max_assoluto_m": 0.001, } + +USCITA = { + # Cartella in cui salvare immagini e animazioni quando il backend + # matplotlib non è interattivo. + "cartella": "output", +} diff --git a/fem.py b/fem.py index 6fefad1..d651e26 100644 --- a/fem.py +++ b/fem.py @@ -1,19 +1,17 @@ # Solutore termico transitorio a elementi finiti sulla mesh shell della # fascetta (vedi mesh.py per la geometria). # -# Differenza rispetto al modello ai volumi finiti di simulate.py: -# - lì il dominio è la sezione x-z (spessore risolto, circonferenza y -# collassata in un'attenuazione gaussiana del flusso e in un sink lineare -# di tipo aletta); -# - qui 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 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 skin depth non entra nel modello shell: conta solo il flusso -# totale assorbito per unità di superficie. +# 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): # @@ -31,21 +29,18 @@ # 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 FEM +from config import ARIA, FASCETTA, FEM, RANDOMIZZAZIONE, SENSORE, SORGENTE from materials import MATERIALI from mesh import genera_mesh -# Helper della cinematica delle sorgenti condivisi con il modello ai volumi -# finiti: la corsa del gruppo di induttori deve essere identica nei due modelli. -from simulate import ( - _intervallo_attivo, - _x_riferimento_finale_m, - _x_riferimento_iniziale_m, -) + +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). @@ -75,6 +70,122 @@ RIGIDEZZA_S = np.array( ) / 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 @@ -124,6 +235,10 @@ def prepara_stato_fem( 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"] @@ -201,6 +316,8 @@ def prepara_stato_fem( "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, } @@ -265,16 +382,29 @@ def campo_iniziale(stato: dict) -> np.ndarray: def simula_campo_fem( - cfg_run: dict, - durata_s: float, 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.""" - fascetta = cfg_run["fascetta"] - aria = cfg_run["aria"] - sorgente = cfg_run["sorgente"] - sensore = cfg_run["sensore"] + """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"] @@ -296,11 +426,14 @@ def simula_campo_fem( 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(T_sensore) + T_lette.append(letta) prossimo_frame_t += dt_frame_s t += dt @@ -315,4 +448,5 @@ def simula_campo_fem( "T_lette": np.array(T_lette), "sorgente": sorgente, "T_ambiente_C": stato["T_ambiente_C"], + "skin_depth_m": stato["skin_depth_m"], } diff --git a/main.py b/main.py index e8b9f2d..3dc1f0d 100644 --- a/main.py +++ b/main.py @@ -2,10 +2,6 @@ # azione eseguire. # # python main.py mesh # disegna la sola mesh a elementi shell -# python main.py simula # genera il dataset di run -# python main.py grafico # grafici del primo run generato -# python main.py anima # animazione 2D della sezione -# python main.py anima3d # animazione 3D isometrica del barattolo # python main.py fem # animazione 3D del campo FEM sulla mesh shell # # Senza argomenti stampa l'elenco delle azioni disponibili. @@ -20,30 +16,6 @@ def _azione_mesh() -> None: plot_mesh.main() -def _azione_simula() -> None: - import simulate - - simulate.main() - - -def _azione_grafico() -> None: - import plot_csv - - plot_csv.main() - - -def _azione_anima() -> None: - import plot_animazione - - plot_animazione.main() - - -def _azione_anima3d() -> None: - import plot_animazione_3d - - plot_animazione_3d.main() - - def _azione_fem() -> None: import plot_animazione_fem @@ -53,17 +25,13 @@ def _azione_fem() -> None: # Chiave da riga di comando -> (funzione, descrizione mostrata nell'help). AZIONI = { "mesh": (_azione_mesh, "Disegna la mesh a elementi shell della fascetta"), - "simula": (_azione_simula, "Genera il dataset di run nella cartella di output"), - "grafico": (_azione_grafico, "Grafici temperatura e flusso del primo run"), - "anima": (_azione_anima, "Animazione 2D del campo di temperatura nella sezione"), - "anima3d": (_azione_anima3d, "Animazione 3D isometrica del barattolo"), "fem": (_azione_fem, "Animazione 3D del campo FEM sulla mesh a elementi shell"), } def main() -> None: parser = argparse.ArgumentParser( - description="Simulatore termico della fascetta riscaldata a induzione.", + description="Analisi termica FEM della fascetta riscaldata a induzione.", formatter_class=argparse.RawTextHelpFormatter, ) parser.add_argument( diff --git a/materials.py b/materials.py index fee5f79..d4806e3 100644 --- a/materials.py +++ b/materials.py @@ -1,4 +1,4 @@ -# Database dei materiali per il simulatore termico. +# Database dei materiali per l'analisi termica. # # Tutte le unità sono SI: # - conducibilita_termica_W_mK diff --git a/plot_animazione.py b/plot_animazione.py deleted file mode 100644 index 8b63413..0000000 --- a/plot_animazione.py +++ /dev/null @@ -1,230 +0,0 @@ -# Animazione della sezione della fascetta durante il passaggio delle sorgenti. -# -# Riproduce la fisica di run_0001 (stesso seed di simulate.py) e mostra: -# - il profilo di flusso termico q(x) sul lato esterno e le sorgenti in moto; -# - il campo di temperatura T(x, z) nella sezione lunghezza × spessore; -# - il sensore infrarosso e la temperatura nel punto osservato. - -import random -from pathlib import Path - -import matplotlib -import matplotlib.pyplot as plt -import numpy as np -from matplotlib.animation import FuncAnimation, PillowWriter - -from config import SIMULAZIONE -from simulate import ( - configurazione_randomizzata, - passo_implicito, - prepara_stato_termico, - profilo_flusso_incidente_W_m2, -) - -# Istante di inizio dei fotogrammi mostrati (la simulazione parte comunque da 0). -T_INIZIO_ANIMAZIONE_S = 0.40 - -# Istante di fine dell'animazione. -T_FINE_ANIMAZIONE_S = 30 - -# Tempo simulato tra un fotogramma e il successivo. -DT_FRAME_S = 0.05 - -# Millisecondi tra i fotogrammi in riproduzione. -INTERVALLO_RIPRODUZIONE_MS = 30 - - -def simula_campi(cfg_run: dict) -> dict: - # Esegue la simulazione fino a T_FINE_ANIMAZIONE_S salvando, a ogni - # fotogramma, campo di temperatura, profilo di flusso e stato del sensore. - fascetta = cfg_run["fascetta"] - aria = cfg_run["aria"] - sorgente = cfg_run["sorgente"] - sensore = cfg_run["sensore"] - - 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 / stato["dx_m"]))) - - 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) - - tempi, campi, flussi, x_riferimenti = [], [], [], [] - T_vere, T_lette = [], [] - - prossimo_frame_t = 0.0 - 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) - T = passo_implicito(stato, T, q_x) - - T_sensore += (T[i_sensore, -1] - T_sensore) * dt / tau_sensore - - if t + 1e-12 >= prossimo_frame_t: - tempi.append(t) - campi.append(T.copy()) - flussi.append(q_x.copy()) - x_riferimenti.append(x_rif) - T_vere.append(T[i_sensore, -1]) - T_lette.append(T_sensore) - prossimo_frame_t += DT_FRAME_S - - t += dt - - return { - "tempi": np.array(tempi), - "campi": campi, - "flussi": flussi, - "x_riferimenti": np.array(x_riferimenti), - "T_vere": np.array(T_vere), - "T_lette": np.array(T_lette), - "x_centri_mm": x_centri * 1000.0, - "spessore_mm": fascetta["spessore_mm"], - "lunghezza_mm": fascetta["lunghezza_mm"], - "x_sensore_mm": sensore["x_mm"], - "sorgente": cfg_run["sorgente"], - } - - -def main() -> None: - rng = random.Random(SIMULAZIONE["seed"]) - cfg_run = configurazione_randomizzata(1, rng) - dati = simula_campi(cfg_run) - - tempi = dati["tempi"] - indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S)) - n_frame = len(tempi) - indice_inizio - - lunghezza_mm = dati["lunghezza_mm"] - spessore_mm = dati["spessore_mm"] - x_vista_mm = (-10.0, lunghezza_mm + 10.0) - - q_max_MW = max(q.max() for q in dati["flussi"]) / 1e6 - T_max = max(c.max() for c in dati["campi"]) - - sorgente = dati["sorgente"] - numero_sorgenti = sorgente.get("numero_sorgenti", 1) - distanza_mm = sorgente.get("distanza_sorgenti_m", 0.0) * 1000.0 - - fig, (ax_flusso, ax_sezione, ax_storia) = plt.subplots( - 3, 1, figsize=(10, 8), height_ratios=[1.0, 1.6, 1.2], - gridspec_kw={"hspace": 0.45}, - ) - fig.suptitle("Sezione della fascetta: sorgenti in transito e sensore") - - # Pannello 1: profilo di flusso sul lato esterno e posizioni delle sorgenti. - linea_flusso, = ax_flusso.plot([], [], color="tab:red") - marker_sorgenti, = ax_flusso.plot( - [], [], "v", color="tab:red", markersize=10, clip_on=False - ) - ax_flusso.annotate( - "verso di marcia", - xy=(0.28, 0.85), xytext=(0.55, 0.85), xycoords="axes fraction", - arrowprops={"arrowstyle": "->", "color": "gray"}, - color="gray", va="center", - ) - ax_flusso.set_xlim(*x_vista_mm) - ax_flusso.set_ylim(0.0, q_max_MW * 1.25) - ax_flusso.set_ylabel("q(x) [MW/m²]") - ax_flusso.grid(True, alpha=0.3) - ax_flusso.set_xticklabels([]) - - # Pannello 2: campo di temperatura nella sezione (z verso il basso, - # origine nel vertice in alto a sinistra come nel modello). - immagine = ax_sezione.imshow( - dati["campi"][indice_inizio].T, - extent=(0.0, lunghezza_mm, spessore_mm, 0.0), - aspect="auto", - cmap="inferno", - vmin=cfg_run["aria"]["temperatura_ambiente_C"], - vmax=T_max, - interpolation="bilinear", - ) - ax_sezione.set_xlim(*x_vista_mm) - ax_sezione.set_ylim(3.2 * spessore_mm, -0.6 * spessore_mm) - ax_sezione.set_ylabel("z [mm]") - ax_sezione.set_xlabel("x [mm]") - # Sensore infrarosso sotto la parete interna (posizione schematica, - # non in scala) con linea di vista tratteggiata. - x_sens = dati["x_sensore_mm"] - ax_sezione.plot([x_sens], [2.4 * spessore_mm], "^", color="tab:blue", markersize=12) - ax_sezione.plot( - [x_sens, x_sens], [1.1 * spessore_mm, 2.1 * spessore_mm], - linestyle="--", color="tab:blue", linewidth=1, - ) - ax_sezione.text( - x_sens + 3, 2.4 * spessore_mm, "sensore IR", - color="tab:blue", va="center", - ) - # La colorbar è agganciata a tutti i pannelli per non restringere solo - # quello della sezione, mantenendo allineati gli assi x. - barra = fig.colorbar( - immagine, ax=(ax_flusso, ax_sezione, ax_storia), pad=0.02, aspect=35 - ) - barra.set_label("T [°C]") - - # Pannello 3: temperatura nel punto osservato dal sensore. - linea_vera, = ax_storia.plot([], [], label="T vera lato interno") - linea_letta, = ax_storia.plot([], [], label="T sensore (con inerzia)") - cursore = ax_storia.axvline(tempi[indice_inizio], color="gray", linewidth=0.8) - ax_storia.set_xlim(0.0, T_FINE_ANIMAZIONE_S) - ax_storia.set_ylim(15.0, max(dati["T_vere"].max(), dati["T_lette"].max()) * 1.08) - ax_storia.set_xlabel("Tempo [s]") - ax_storia.set_ylabel("T [°C]") - ax_storia.legend(loc="upper left") - ax_storia.grid(True, alpha=0.3) - - testo_tempo = ax_flusso.set_title(f"t = {tempi[indice_inizio]:.3f} s", loc="right") - - def aggiorna(frame: int): - k = indice_inizio + frame - t = tempi[k] - - linea_flusso.set_data(dati["x_centri_mm"], dati["flussi"][k] / 1e6) - - x_sorgenti_mm = ( - dati["x_riferimenti"][k] * 1000.0 - + np.arange(numero_sorgenti) * distanza_mm - ) - visibili = (x_sorgenti_mm >= x_vista_mm[0]) & (x_sorgenti_mm <= x_vista_mm[1]) - marker_sorgenti.set_data( - x_sorgenti_mm[visibili], - np.full(int(visibili.sum()), q_max_MW * 1.12), - ) - - immagine.set_data(dati["campi"][k].T) - - linea_vera.set_data(tempi[: k + 1], dati["T_vere"][: k + 1]) - linea_letta.set_data(tempi[: k + 1], dati["T_lette"][: k + 1]) - cursore.set_xdata([t, t]) - testo_tempo.set_text(f"t = {t:.3f} s") - - return ( - linea_flusso, marker_sorgenti, immagine, - linea_vera, linea_letta, cursore, testo_tempo, - ) - - animazione = FuncAnimation( - fig, aggiorna, frames=n_frame, interval=INTERVALLO_RIPRODUZIONE_MS, blit=False - ) - - # Se il backend non è interattivo si salva una GIF invece di mostrare la finestra. - if matplotlib.get_backend().lower() == "agg": - percorso = Path("dataset") / "animazione_sezione.gif" - animazione.save( - percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS) - ) - print(f"Backend non interattivo: animazione salvata in {percorso}") - return - - plt.show() - - -if __name__ == "__main__": - main() diff --git a/plot_animazione_3d.py b/plot_animazione_3d.py deleted file mode 100644 index 9207912..0000000 --- a/plot_animazione_3d.py +++ /dev/null @@ -1,150 +0,0 @@ -# Animazione 3D del barattolo in vista isometrica: colormap della temperatura -# sulla superficie esterna. -# -# Il modello risolve solo la sezione x-z (vedi CLAUDE.md): la coordinata -# circonferenziale y non è simulata, è collassata in un'attenuazione -# gaussiana del flusso. Per estrudere il campo attorno alla circonferenza si -# riusa la stessa gaussiana: la sovratemperatura rispetto al piano y=0 (dove -# si trova il sensore) viene scalata in funzione della distanza -# circonferenziale dal percorso delle sorgenti. È quindi una ricostruzione -# visiva, non un calcolo di diffusione in y. - -import random - -import matplotlib -import matplotlib.pyplot as plt -import numpy as np -from matplotlib import cm -from matplotlib.animation import FuncAnimation, PillowWriter -from pathlib import Path - -from config import FASCETTA, SIMULAZIONE -from plot_animazione import DT_FRAME_S as DT_FRAME_S_BASE -from plot_animazione import T_FINE_ANIMAZIONE_S, T_INIZIO_ANIMAZIONE_S, simula_campi -from simulate import configurazione_randomizzata - -# Tempo simulato tra un fotogramma e il successivo dell'animazione 3D. Più -# largo di DT_FRAME_S_BASE perché ricreare una superficie 3D a ogni -# fotogramma è più costoso della semplice imshow 2D. -DT_FRAME_S = 0.15 - -# Millisecondi tra i fotogrammi in riproduzione. -INTERVALLO_RIPRODUZIONE_MS = 60 - -# Numero di punti lungo la circonferenza per l'estrusione della superficie. -N_THETA = 72 - - -def attenuazione_circonferenziale( - y_m: np.ndarray, offset_y_m: float, sigma_m: float -) -> np.ndarray: - # Fattore che scala la sovratemperatura (T - T_ambiente) del piano y=0 - # in funzione della distanza circonferenziale y dal piano stesso, - # normalizzato in modo da valere 1 in y=0. - sigma = max(sigma_m, 1e-9) - esponente = -0.5 * ((y_m - offset_y_m) ** 2 - offset_y_m ** 2) / (sigma * sigma) - return np.exp(esponente) - - -def main() -> None: - rng = random.Random(SIMULAZIONE["seed"]) - cfg_run = configurazione_randomizzata(1, rng) - dati = simula_campi(cfg_run) - - tempi = dati["tempi"] - indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S)) - passo = max(1, round(DT_FRAME_S / DT_FRAME_S_BASE)) - indici_frame = list(range(indice_inizio, len(tempi), passo)) - - T_ambiente = cfg_run["aria"]["temperatura_ambiente_C"] - sorgente = cfg_run["sorgente"] - sigma_m = sorgente["sigma_punto_m"] - offset_y_m = sorgente["offset_y_percorso_m"] - numero_sorgenti = sorgente.get("numero_sorgenti", 1) - distanza_m = sorgente.get("distanza_sorgenti_m", 0.0) - - raggio_m = (FASCETTA["diametro_mm"] / 1000.0) / 2.0 - x_centri_m = dati["x_centri_mm"] / 1000.0 - lunghezza_m = dati["lunghezza_mm"] / 1000.0 - - theta = np.linspace(-np.pi, np.pi, N_THETA) - y_circ_m = theta * raggio_m - attenuazione = attenuazione_circonferenziale(y_circ_m, offset_y_m, sigma_m) - - Xm, Thetam = np.meshgrid(x_centri_m, theta) - Ym = raggio_m * np.sin(Thetam) - Zm = raggio_m * np.cos(Thetam) - - T_max = max(c[:, 0].max() for c in dati["campi"]) - # vmin più basso della temperatura ambiente reale: altrimenti le zone - # fredde cadrebbero sul nero puro di "inferno" e, essendo lo shading - # moltiplicativo, nessuna illuminazione basterebbe a renderle visibili. - norm = matplotlib.colors.Normalize( - vmin=T_ambiente - 0.4 * (T_max - T_ambiente), vmax=T_max - ) - cmap = matplotlib.colormaps["inferno"] - lightsource = matplotlib.colors.LightSource(azdeg=315, altdeg=45) - - fig = plt.figure(figsize=(9, 7)) - ax = fig.add_subplot(projection="3d") - ax.view_init(elev=35.264, azim=45) - ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m)) - ax.set_xlabel("x [m]") - ax.set_axis_off() - - # La colorbar mostra il range reale delle temperature: il norm esteso - # verso il basso serve solo a schiarire il colore di base della - # superficie fredda, non deve comparire nella scala mostrata all'utente. - norm_colorbar = matplotlib.colors.Normalize(vmin=T_ambiente, vmax=T_max) - mappabile = cm.ScalarMappable(cmap=cmap, norm=norm_colorbar) - mappabile.set_array([]) - barra = fig.colorbar(mappabile, ax=ax, shrink=0.6, pad=0.05) - barra.set_label("T [°C]") - - def disegna_frame(k: int): - ax.cla() - ax.view_init(elev=35.264, azim=45) - ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m)) - ax.set_axis_off() - - T_lato_esterno = dati["campi"][k][:, 0] - T_superficie = T_ambiente + (T_lato_esterno[None, :] - T_ambiente) * attenuazione[:, None] - colori = cmap(norm(T_superficie)) - - ax.plot_surface( - Xm, Ym, Zm, facecolors=colori, rstride=1, cstride=1, - antialiased=False, shade=True, lightsource=lightsource, linewidth=0, - ) - - x_sorgenti_m = dati["x_riferimenti"][k] + np.arange(numero_sorgenti) * distanza_m - visibili = (x_sorgenti_m >= 0.0) & (x_sorgenti_m <= lunghezza_m) - if visibili.any(): - theta_sorgente = offset_y_m / raggio_m - ax.scatter( - x_sorgenti_m[visibili], - np.full(int(visibili.sum()), raggio_m * np.sin(theta_sorgente) * 1.05), - np.full(int(visibili.sum()), raggio_m * np.cos(theta_sorgente) * 1.05), - color="cyan", s=25, depthshade=False, - ) - - ax.set_title(f"t = {tempi[k]:.3f} s") - return () - - animazione = FuncAnimation( - fig, disegna_frame, frames=indici_frame, - interval=INTERVALLO_RIPRODUZIONE_MS, blit=False, - ) - - if matplotlib.get_backend().lower() == "agg": - percorso = Path("dataset") / "animazione_3d.gif" - animazione.save( - percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS) - ) - print(f"Backend non interattivo: animazione salvata in {percorso}") - return - - plt.show() - - -if __name__ == "__main__": - main() diff --git a/plot_animazione_fem.py b/plot_animazione_fem.py index 27af358..806cafe 100644 --- a/plot_animazione_fem.py +++ b/plot_animazione_fem.py @@ -1,17 +1,15 @@ -# Animazione 3D isometrica del barattolo con il campo di temperatura calcolato +# Animazione 3D isometrica della fascetta con il campo di temperatura calcolato # a elementi finiti sulla mesh shell (fem.py). # -# È l'equivalente di plot_animazione_3d.py ma senza ricostruzioni: lì il campo -# circonferenziale era estruso riusando la gaussiana della sorgente, qui la -# coordinata circonferenziale è una direzione risolta del modello FEM, quindi -# la temperatura dipinta su ogni elemento è quella effettivamente calcolata. +# La coordinata circonferenziale è una direzione risolta del modello FEM, +# quindi la temperatura dipinta su ogni elemento è quella effettivamente +# calcolata, senza ricostruzioni. # # Ogni elemento shell è disegnato come faccia piana con colore pari alla media # dei suoi quattro valori nodali. L'ombreggiatura è calcolata a mano dalla # normale radiale dell'elemento: Poly3DCollection ombreggia solo alla # creazione, mentre qui i colori delle facce cambiano a ogni fotogramma. -import random from pathlib import Path import matplotlib @@ -21,11 +19,13 @@ from matplotlib import cm from matplotlib.animation import FuncAnimation, PillowWriter from mpl_toolkits.mplot3d.art3d import Poly3DCollection -from config import SIMULAZIONE +from config import USCITA from fem import simula_campo_fem from mesh import genera_mesh, riepilogo_mesh -from plot_animazione import T_FINE_ANIMAZIONE_S, T_INIZIO_ANIMAZIONE_S -from simulate import configurazione_randomizzata + +# Istante del primo fotogramma mostrato: i primi istanti sono ancora uniformi +# alla temperatura ambiente. +T_INIZIO_ANIMAZIONE_S = 0.40 # Tempo simulato tra un fotogramma e il successivo. DT_FRAME_S = 0.15 @@ -65,16 +65,8 @@ def main() -> None: mesh_dati = genera_mesh() print(riepilogo_mesh(mesh_dati)) - rng = random.Random(SIMULAZIONE["seed"]) - cfg_run = configurazione_randomizzata(1, rng) - print("Integrazione FEM in corso...") - dati = simula_campo_fem( - cfg_run, - durata_s=T_FINE_ANIMAZIONE_S, - dt_frame_s=DT_FRAME_S, - mesh_dati=mesh_dati, - ) + dati = simula_campo_fem(dt_frame_s=DT_FRAME_S, mesh_dati=mesh_dati) tempi = dati["tempi"] indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S)) @@ -84,7 +76,8 @@ def main() -> None: T_max = max(campo.max() for campo in dati["campi"]) print( f"T massima nodale: {T_max:.1f} °C — " - f"T massima nel punto del sensore: {dati['T_vere'].max():.1f} °C" + f"T massima nel punto del sensore: {dati['T_vere'].max():.1f} °C\n" + f"Skin depth (diagnostica): {1000.0 * dati['skin_depth_m']:.3f} mm" ) elementi = mesh_dati["elementi"] @@ -133,17 +126,17 @@ def main() -> None: # Pannello inferiore: storia della temperatura nel punto osservato dal # sensore. T_vere è il valore nodale calcolato dal FEM, T_lette lo stesso - # segnale filtrato dall'inerzia del primo ordine del sensore. - # Con la costante di tempo del sensore le due curve sono quasi - # sovrapposte: la prima è tracciata spessa e trasparente perché la seconda, - # tratteggiata, resti leggibile sopra di essa. + # segnale degradato dal sensore reale (inerzia, rumore, quantizzazione). + # A parte il rumore le due curve sono quasi sovrapposte: la prima è + # tracciata spessa e trasparente perché la seconda resti leggibile sopra + # di essa. linea_vera, = ax_storia.plot( [], [], color="tab:blue", linewidth=3.0, alpha=0.4, label="T nodo del sensore (FEM)", ) linea_letta, = ax_storia.plot( - [], [], color="tab:red", linewidth=1.2, linestyle="--", - label="T misurata dal sensore (con inerzia)", + [], [], color="tab:red", linewidth=1.0, + label="T misurata dal sensore (inerzia + rumore + quantizzazione)", ) cursore = ax_storia.axvline(tempi[indice_inizio], color="gray", linewidth=0.8) ax_storia.set_xlim(0.0, tempi[-1]) @@ -189,7 +182,7 @@ def main() -> None: ) if matplotlib.get_backend().lower() == "agg": - cartella = Path("dataset") + cartella = Path(USCITA["cartella"]) cartella.mkdir(parents=True, exist_ok=True) percorso = cartella / "animazione_fem.gif" animazione.save( diff --git a/plot_csv.py b/plot_csv.py deleted file mode 100644 index 287a869..0000000 --- a/plot_csv.py +++ /dev/null @@ -1,48 +0,0 @@ -from pathlib import Path - -import matplotlib -import matplotlib.pyplot as plt -import pandas as pd - - -def main() -> None: - percorso_csv = Path("dataset/run_0001.csv") - if not percorso_csv.exists(): - raise FileNotFoundError( - "File CSV non trovato. Esegui prima `python simulate.py`." - ) - - df = pd.read_csv(percorso_csv) - - fig_temperature = plt.figure() - plt.plot(df["tempo_s"], df["T_vera_lato_sensore_C"], label="Temperatura vera lato interno") - plt.plot(df["tempo_s"], df["T_misurata_sensore_C"], label="Temperatura misurata dal sensore") - plt.xlabel("Tempo [s]") - plt.ylabel("Temperatura [°C]") - plt.title("Sensore fisso, sorgenti in moto") - plt.legend() - plt.grid(True) - plt.tight_layout() - - fig_flusso = plt.figure() - plt.plot(df["tempo_s"], df["flusso_termico_sorgente_W_m2"]) - plt.xlabel("Tempo [s]") - plt.ylabel("Flusso termico efficace [W/m²]") - plt.title("Flusso termico nel punto x osservato dal sensore") - plt.grid(True) - plt.tight_layout() - - # Se il backend non è interattivo (nessun display o backend GUI - # disponibile), plt.show() non aprirebbe nulla: si salvano i PNG. - if matplotlib.get_backend().lower() == "agg": - for fig, nome in [(fig_temperature, "temperature"), (fig_flusso, "flusso")]: - percorso = Path("dataset") / f"grafico_{nome}.png" - fig.savefig(percorso, dpi=150) - print(f"Backend non interattivo: grafico salvato in {percorso}") - return - - plt.show() - - -if __name__ == "__main__": - main() diff --git a/plot_mesh.py b/plot_mesh.py index 6c9201c..6ccfae8 100644 --- a/plot_mesh.py +++ b/plot_mesh.py @@ -7,6 +7,7 @@ import matplotlib import matplotlib.pyplot as plt from mpl_toolkits.mplot3d.art3d import Poly3DCollection +from config import USCITA from mesh import genera_mesh, riepilogo_mesh @@ -45,7 +46,7 @@ def main() -> None: ) if matplotlib.get_backend().lower() == "agg": - cartella = Path("dataset") + cartella = Path(USCITA["cartella"]) cartella.mkdir(parents=True, exist_ok=True) percorso = cartella / "mesh.png" fig.savefig(percorso, dpi=150) diff --git a/simulate.py b/simulate.py deleted file mode 100644 index 7b482b0..0000000 --- a/simulate.py +++ /dev/null @@ -1,486 +0,0 @@ -import csv -import math -import os -import random -import shutil -from concurrent.futures import ProcessPoolExecutor -from copy import deepcopy -from pathlib import Path - -import numpy as np -import scipy.sparse as sp -from scipy.sparse.linalg import splu - -from config import ARIA, FASCETTA, 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, 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 profilo_flusso_incidente_W_m2( - sorgente: dict, - x_sensore_m: float, - t_s: float, - x_centri_m: np.ndarray, -) -> tuple[float, np.ndarray]: - # Restituisce x_sorgente_m (posizione della sorgente di riferimento) e - # il profilo di flusso termico efficace q(x) [W/m²] sul lato esterno, - # somma dei contributi di tutte le sorgenti attive. - # - # Ogni sorgente in moto ha un'impronta gaussiana lungo x, valutata sui - # centri cella della sezione. L'offset circonferenziale y non è risolto - # spazialmente: entra come attenuazione gaussiana del flusso. - 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) - dy = sorgente["offset_y_percorso_m"] - sigma = sorgente["sigma_punto_m"] - - q_picco = sorgente["flusso_termico_picco_W_m2"] * sorgente["efficienza_riscaldamento"] - attenuazione_y = math.exp(-0.5 * (dy * dy) / (sigma * sigma)) - - q_x = np.zeros_like(x_centri_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 = x_centri_m - x_i - q_x += q_picco * attenuazione_y * np.exp(-0.5 * (dx * dx) / (sigma * sigma)) - - return x_riferimento, q_x - - -def profilo_deposizione_z_1_m( - z_centri_m: np.ndarray, - spessore_m: float, - skin_depth_m: float, -) -> np.ndarray: - # Profilo di deposizione volumetrica del flusso superficiale [1/m]: - # - # p(z) = exp(-z / delta) / (delta * (1 - exp(-spessore / delta))) - # - # Normalizzato in modo che integrale_0^spessore p(z) dz = 1, così che - # q_vol(x, z) = q(x) * p(z) conservi il flusso superficiale. - delta = max(skin_depth_m, 1e-9) - normalizzazione = delta * (1.0 - math.exp(-spessore_m / delta)) - return np.exp(-z_centri_m / delta) / normalizzazione - - -def _laplaciano_1d(n: int) -> sp.spmatrix: - # Operatore alle differenze -T'' su n celle con bordi adiabatici (Neumann). - diagonale = np.full(n, 2.0) - diagonale[0] = 1.0 - diagonale[-1] = 1.0 - fuori = -np.ones(n - 1) - return sp.diags([fuori, diagonale, fuori], [-1, 0, 1]) - - -def costruisci_solutore_implicito_2d( - n_x: int, - n_z: int, - dt_s: float, - dx_m: float, - dz_m: float, - spessore_m: float, - materiale: dict, - h_esterno_W_m2K: float, - h_interno_W_m2K: float, - h_bordi_W_m2K: 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). 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. Il termine è ricavato - # considerando l'intero volume della fascia cilindrica (equazione - # dell'aletta): il calore conduce lungo y attraverso l'intero spessore - # mentre le superfici esterna e interna dell'intero cilindro scambiano - # per convezione, dando un sink distribuito uniformemente su ogni cella - # q_y = -(h_esterno + h_interno) / spessore * (T - T_amb). - k = materiale["conducibilita_termica_W_mK"] - rho = materiale["densita_kg_m3"] - cp = materiale["calore_specifico_J_kgK"] - alpha = k / (rho * cp) - - r_x = alpha * dt_s / (dx_m * dx_m) - 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 = (h_esterno_W_m2K + h_interno_W_m2K) * dt_s / (rho * cp * spessore_m) - - n = n_x * n_z - 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(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 - - 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": (aria["h_esterno_W_m2K"] + aria["h_interno_W_m2K"]) * dt / (rho * cp * spessore), - "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, - spessore_m=spessore, - 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"], - ), - } - - -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 - return round(valore / passo) * passo - - -def configurazione_randomizzata(indice_run: int, rng: random.Random) -> dict: - fascetta = deepcopy(FASCETTA) - 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}", - "fascetta": fascetta, - "aria": aria, - "sorgente": sorgente, - "sensore": sensore, - } - - -def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: - fascetta = cfg_run["fascetta"] - aria = cfg_run["aria"] - sorgente = cfg_run["sorgente"] - sensore = cfg_run["sensore"] - - nome_materiale = fascetta["materiale"] - - 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 / stato["dx_m"]))) - - dt = stato["dt_s"] - durata = SIMULAZIONE["durata_s"] - periodo_campionamento = 1.0 / SIMULAZIONE["frequenza_campionamento_hz"] - - T = np.full((n_x, n_z), aria["temperatura_ambiente_C"], dtype=float) - T_sensore = T[i_sensore, -1] - - prossimo_campione_t = 0.0 - T_vera_max = T[i_sensore, -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_x = profilo_flusso_incidente_W_m2( - sorgente, x_sensore, t, x_centri - ) - T = passo_implicito(stato, T, q_x) - - # Temperatura vera della superficie interna nel punto osservato - # dal sensore infrarosso. - T_vera_lato_sensore = T[i_sensore, -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_x[i_sensore]:.6f}", - f"{skin_depth:.9e}", - f"{T_vera_lato_sensore:.6f}", - f"{misurata:.6f}", - f"{T[i_sensore, 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, - "diametro_m": fascetta["diametro_mm"] / 1000.0, - "lunghezza_m": fascetta["lunghezza_mm"] / 1000.0, - "spessore_m": fascetta["spessore_mm"] / 1000.0, - "n_nodi_x": n_x, - "n_nodi_z": n_z, - "durata_s": durata, - "frequenza_campionamento_hz": SIMULAZIONE["frequenza_campionamento_hz"], - "dt_interno_s": dt, - "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, - "distanza_sensore_parete_m": sensore["distanza_parete_mm"] / 1000.0, - "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 _esegui_run(indice_e_seme: tuple[int, int]) -> dict: - # Ogni run riceve un seme indipendente derivato dal seed globale, così - # l'esecuzione in parallelo resta riproducibile indipendentemente - # dall'ordine in cui i processi la completano. - indice, seme = indice_e_seme - rng = random.Random(seme) - cfg_run = configurazione_randomizzata(indice, rng) - cartella_output = Path(SIMULAZIONE["cartella_output"]) - percorso_csv = cartella_output / f"{cfg_run['id_run']}.csv" - return simula_singolo(cfg_run, percorso_csv, rng) - - -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_semi = random.Random(SIMULAZIONE["seed"]) - num_run = SIMULAZIONE["num_run"] - semi = [rng_semi.randrange(2**63) for _ in range(num_run)] - - num_processi = SIMULAZIONE["num_processi"] or os.cpu_count() or 1 - with ProcessPoolExecutor(max_workers=num_processi) as executor: - righe_metadata = list( - executor.map(_esegui_run, enumerate(semi, start=1)) - ) - - 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()