8 Commits
Author SHA1 Message Date
davideandClaude Opus 5 97ebab844e Esporta in CSV la temperatura simulata e la lettura del sensore
Nuova azione `python main.py csv`: esporta_csv.py integra il campo FEM e scrive
output/csv/run_XXXX.csv con la temperatura vera della parete nel punto
osservato, la lettura del sensore reale e il massimo nodale istante per
istante, più un metadata.csv con una riga per analisi. La cartella è ricreata a
ogni esecuzione.

Le righe sono campionate a FEM["frequenza_campionamento_hz"],
indipendente dal passo di integrazione, e FEM["num_run"] analisi vengono
esportate di seguito: con la randomizzazione abilitata ognuna ha parametri
diversi ma riproducibili dallo stesso seed.

simula_campo_fem accetta salva_campi=False, perché all'export servono solo le
serie scalari e non i campi nodali di ogni campione, e restituisce anche la
temperatura nodale massima per campione.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-03 11:44:58 +02:00
davideandClaude Opus 5 1631d6c94e Rimuove il modello ai volumi finiti lasciando la sola analisi FEM
Il progetto conteneva due discretizzazioni indipendenti della stessa fisica:
i volumi finiti sulla sezione x-z di simulate.py (con generazione del dataset
CSV e le sue animazioni 2D/3D) e gli elementi finiti sulla mesh shell di
fem.py. Resta solo il secondo, che risolve la circonferenza invece di
collassarla in un'attenuazione gaussiana.

Eliminati simulate.py, plot_csv.py, plot_animazione.py e plot_animazione_3d.py.
Il FEM dipendeva da simulate.py per la cinematica del gruppo di sorgenti e per
la configurazione del run, quindi fem.py assorbe _x_riferimento_iniziale_m,
_x_riferimento_finale_m, _intervallo_attivo e configurazione_randomizzata; la
lettura del sensore (inerzia, rumore, quantizzazione) è ora applicata dentro
simula_campo_fem e la skin depth resta calcolata dalla frequenza dell'induttore
come sola grandezza diagnostica, dato che nel modello shell lo spessore è
collassato.

config.py perde i parametri che servivano solo ai volumi finiti (griglia
n_nodi_x/n_nodi_z, campionamento e cartella del dataset, numero di processi) e
guadagna FEM["durata_s"], FEM["seed"] e USCITA["cartella"]: le immagini e le
animazioni vanno in output/ invece che in dataset/.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-03 11:44:58 +02:00
davideandClaude Opus 5 6a55427b60 Valuta il flusso nodale sfruttando la separabilita della gaussiana
L'impronta della sorgente e separabile in (x, s) e il fattore
circonferenziale e costante durante il run: viene precalcolato una volta in
prepara_stato_fem, mentre il fattore assiale e valutato sui soli n_x + 1 nodi
distinti in x ed espanso con un prodotto esterno. Il risultato e identico
(differenza max 1e-104) e il costo del termine noto scende da 108 a 23 us per
passo, cioe -31% sul tempo totale di integrazione.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-03 10:08:16 +02:00
davideandClaude Opus 5 840f6fd09e Aggiunge il pannello T-t del sensore all'animazione FEM
Sotto la vista 3D compare la storia della temperatura nel nodo osservato dal
sensore, con il cursore sull'istante corrente: la curva nodale calcolata dal
FEM e la stessa filtrata dall'inerzia del sensore.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-03 09:43:17 +02:00
davideandClaude Opus 5 0f46f3d479 Aggiunge il solutore termico FEM sulla mesh shell e la sua animazione 3D
fem.py risolve l'equazione del calore transitoria sulla superficie media del
cilindro con elementi shell quadrangolari bilineari: la circonferenza diventa
una direzione risolta (niente più attenuazione gaussiana né sink di aletta),
mentre lo spessore è collassato perché con parete da 0.18 mm il numero di
Biot è ~1e-7. Eulero implicito con matrice costante fattorizzata una volta
per run; le matrici di elemento sono identiche per tutti gli elementi e
l'assemblaggio è vettorizzato.

plot_animazione_fem.py dipinge il campo calcolato direttamente sugli elementi
della mesh, con ombreggiatura ricavata dalla normale radiale perché i colori
delle facce cambiano a ogni fotogramma.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-03 09:29:35 +02:00
davideandClaude Opus 5 14226234bd Aggiunge la mesh a elementi shell della fascetta e il main.py
Introduce mesh.py, che costruisce una griglia strutturata di quadrilateri
a 4 nodi sulla superficie media del cilindro, chiusa lungo la
circonferenza, come base per la futura analisi FEM. plot_mesh.py la
visualizza in vista isometrica e main.py raccoglie tutte le azioni del
progetto dietro un unico punto di ingresso.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-03 09:10:11 +02:00
davide 560c6a6a62 Aggiunge animazione 3D isometrica del barattolo con colormap di temperatura
Estrude il campo x-z della superficie esterna lungo la circonferenza
tramite la stessa gaussiana usata per l'attenuazione del flusso, con
shading per rendere visibile la curvatura del cilindro anche nelle
zone a temperatura ambiente.
2026-07-06 11:59:34 +02:00
davide 288ae81b3e Parallelizza generazione run su più processi CPU
Ogni run usa un seme RNG indipendente derivato dal seed globale, così
l'esecuzione su ProcessPoolExecutor resta riproducibile.
2026-07-06 11:38:33 +02:00
14 changed files with 1276 additions and 976 deletions
+1 -1
View File
@@ -1,3 +1,3 @@
.venv/
__pycache__/
dataset/
output/
+20 -26
View File
@@ -18,45 +18,39 @@ source .venv/bin/activate
# Installare le dipendenze
pip install -r requirements.txt
# Generare il dataset (scrive dataset/run_XXXX.csv + dataset/metadata.csv)
python simulate.py
# Visualizzare il primo run
python plot_csv.py
# 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 fem # animazione 3D del campo FEM sulla mesh shell
python main.py csv # esporta output/csv/run_XXXX.csv + metadata.csv
```
Ogni modulo resta eseguibile anche direttamente (`python plot_mesh.py`, `python plot_animazione_fem.py`, `python esporta_csv.py`).
Attivare sempre il venv (`source .venv/bin/activate`) prima di eseguire qualsiasi comando Python.
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:**
1. `config.py`tutti i parametri configurabili (dizionari SIMULAZIONE, FASCETTA, 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
1. `main.py`dispatcher da riga di comando che seleziona l'azione (mesh, fem, csv)
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` ed esportato in CSV da `esporta_csv.py`
**Pipeline fisica dentro `simula_singolo()` in [simulate.py](simulate.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 coordinate globali x = asse, y = R·sin(theta), z = R·cos(theta), con theta = 0 sul piano del sensore.
- 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
**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`).
`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.
**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`.
**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.
**Export CSV (`esporta_csv.py`):** `FEM["num_run"]` analisi (una per file `output/csv/run_XXXX.csv`), righe campionate a `FEM["frequenza_campionamento_hz"]` indipendentemente da `FEM["dt_s"]`. Colonne: `id_run, tempo_s, x_sorgente_m, T_vera_sensore_C, T_misurata_sensore_C, T_massima_fascetta_C, T_ambiente_C, offset_y_sorgente_m, velocita_m_s, sigma_punto_m, flusso_picco_W_m2, skin_depth_m, materiale`. `metadata.csv` ha una riga per analisi con tutti i parametri e i picchi. L'export chiama `simula_campo_fem(salva_campi=False)`: usa solo le serie scalari, i campi nodali non vengono accumulati. La cartella `output/csv` è cancellata e ricreata a ogni esecuzione (`shutil.rmtree`).
## Convenzioni su `config.py`
@@ -64,8 +58,8 @@ 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, e la conduzione lungo y è un termine di scambio lineare verso la temperatura ambiente; il diametro è registrato solo come geometria del setup.
- 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).
+156 -145
View File
@@ -1,113 +1,135 @@
# 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 o esportato in CSV come serie temporale (temperatura vera
della parete e lettura del sensore).
## 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
esporta_csv.py export CSV della serie temporale
output/ immagini e GIF salvate quando il backend non è interattivo
output/csv/ CSV esportati (ricreata a ogni esecuzione)
```
## Installazione
@@ -121,38 +143,67 @@ 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
python main.py csv # esporta la serie temporale in output/csv/
```
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`, `python esporta_csv.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"]`.
## Output CSV
`python main.py csv` cancella e ricrea `output/csv/`, poi esporta `FEM["num_run"]` analisi
(una per file `run_XXXX.csv`) più un `metadata.csv` con una riga per analisi. Le righe
sono campionate a `FEM["frequenza_campionamento_hz"]`, indipendente dal passo di
integrazione.
### `output/csv/run_XXXX.csv` — serie temporale
| Colonna | Significato |
|-------------------------|-----------------------------------------------------------------|
| `id_run` | identificativo dell'analisi |
| `tempo_s` | tempo simulato |
| `x_sorgente_m` | posizione della sorgente di riferimento del gruppo |
| `T_vera_sensore_C` | temperatura vera della parete nel punto osservato (valore nodale FEM) |
| `T_misurata_sensore_C` | lettura del sensore (inerzia + rumore + quantizzazione) |
| `T_massima_fascetta_C` | temperatura nodale massima su tutta la fascetta in quell'istante |
| `T_ambiente_C` | temperatura ambiente dell'analisi |
| `offset_y_sorgente_m`, `velocita_m_s`, `sigma_punto_m`, `flusso_picco_W_m2` | parametri effettivi dell'analisi (variano se la randomizzazione è attiva) |
| `skin_depth_m` | skin depth diagnostica |
| `materiale` | chiave del materiale |
### `output/csv/metadata.csv` — una riga per analisi
Tutti i parametri effettivi (geometria, mesh, coefficienti di scambio, sorgenti, sensore)
e i valori di picco: `T_vera_max_sensore_C`, `T_misurata_max_sensore_C`,
`T_massima_fascetta_C`.
## 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, seed, campionamento CSV, numero di analisi |
| `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 +214,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,150,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.
+52 -55
View File
@@ -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,28 +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",
}
FASCETTA = {
# Diametro della fascetta [mm].
"diametro_mm": 70.0,
@@ -65,28 +34,49 @@ 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",
}
MESH = {
# Numero di elementi shell lungo l'asse x (lunghezza della fascetta).
"n_elementi_x": 40,
# Numero di elementi shell lungo la circonferenza.
# La mesh è chiusa su se stessa: non c'è una riga di nodi duplicata.
"n_elementi_circonferenza": 48,
}
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,
# Frequenza di campionamento delle righe del CSV esportato.
# Esempio: 10 Hz significa una riga ogni 0.1 s.
"frequenza_campionamento_hz": 10.0,
# Numero di analisi da esportare in CSV, una per file.
# Con la randomizzazione abilitata ogni analisi ha parametri diversi.
"num_run": 1,
}
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,
}
@@ -123,7 +113,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.
@@ -138,7 +129,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,
@@ -157,7 +148,7 @@ SENSORE = {
}
RANDOMIZZAZIONE = {
# Se abilitata, ogni run varia leggermente alcuni parametri.
# Se abilitata, ogni analisi varia leggermente alcuni parametri.
"abilitata": False,
# Deviazioni standard relative.
@@ -174,3 +165,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",
}
+153
View File
@@ -0,0 +1,153 @@
# Export CSV della serie temporale calcolata dal solutore FEM (fem.py).
#
# Per ogni analisi viene scritto un file run_XXXX.csv con la temperatura vera
# della parete nel punto osservato dal sensore e la lettura del sensore reale
# (inerzia, rumore, quantizzazione), più i parametri effettivi dell'analisi.
# `metadata.csv` riassume un'analisi per riga.
#
# Le righe sono campionate a FEM["frequenza_campionamento_hz"], indipendente
# dal passo di integrazione FEM["dt_s"]. I campi nodali non vengono accumulati:
# servono solo le serie scalari.
import csv
import random
import shutil
from pathlib import Path
from config import FEM, USCITA
from fem import configurazione_randomizzata, simula_campo_fem
from mesh import genera_mesh, riepilogo_mesh
INTESTAZIONE = [
"id_run",
"tempo_s",
"x_sorgente_m",
"T_vera_sensore_C",
"T_misurata_sensore_C",
"T_massima_fascetta_C",
"T_ambiente_C",
"offset_y_sorgente_m",
"velocita_m_s",
"sigma_punto_m",
"flusso_picco_W_m2",
"skin_depth_m",
"materiale",
]
def esporta_run(
id_run: str,
percorso_csv: Path,
mesh_dati: dict,
rng: random.Random,
) -> dict:
"""Integra un'analisi FEM, ne scrive il CSV e restituisce i suoi metadati."""
cfg = configurazione_randomizzata(rng)
dati = simula_campo_fem(
dt_frame_s=1.0 / FEM["frequenza_campionamento_hz"],
cfg=cfg,
mesh_dati=mesh_dati,
rng=rng,
salva_campi=False,
)
fascetta = cfg["fascetta"]
aria = cfg["aria"]
sorgente = cfg["sorgente"]
sensore = cfg["sensore"]
skin_depth = dati["skin_depth_m"]
percorso_csv.parent.mkdir(parents=True, exist_ok=True)
with percorso_csv.open("w", newline="") as f:
writer = csv.writer(f)
writer.writerow(INTESTAZIONE)
for k, t in enumerate(dati["tempi"]):
writer.writerow([
id_run,
f"{t:.6f}",
f"{dati['x_riferimenti'][k]:.9f}",
f"{dati['T_vere'][k]:.6f}",
f"{dati['T_lette'][k]:.6f}",
f"{dati['T_massime'][k]:.6f}",
f"{aria['temperatura_ambiente_C']:.6f}",
f"{sorgente['offset_y_percorso_m']:.9f}",
f"{sorgente['velocita_m_s']:.9f}",
f"{sorgente['sigma_punto_m']:.9f}",
f"{sorgente['flusso_termico_picco_W_m2']:.6f}",
f"{skin_depth:.9e}",
fascetta["materiale"],
])
return {
"id_run": id_run,
"file_csv": percorso_csv.name,
"materiale": fascetta["materiale"],
"diametro_m": fascetta["diametro_mm"] / 1000.0,
"lunghezza_m": fascetta["lunghezza_mm"] / 1000.0,
"spessore_m": fascetta["spessore_mm"] / 1000.0,
"n_elementi_x": mesh_dati["n_elementi_x"],
"n_elementi_circonferenza": mesh_dati["n_elementi_circonferenza"],
"durata_s": FEM["durata_s"],
"dt_s": FEM["dt_s"],
"frequenza_campionamento_hz": FEM["frequenza_campionamento_hz"],
"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": dati["stato"]["x_sensore_m"],
"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_sensore_C": float(dati["T_vere"].max()),
"T_misurata_max_sensore_C": float(dati["T_lette"].max()),
"T_massima_fascetta_C": float(dati["T_massime"].max()),
}
def main() -> None:
mesh_dati = genera_mesh()
print(riepilogo_mesh(mesh_dati))
# La cartella dei CSV è ricreata da zero a ogni esecuzione, così non
# restano run di esecuzioni precedenti con parametri diversi.
cartella = Path(USCITA["cartella"]) / "csv"
if cartella.exists():
shutil.rmtree(cartella)
cartella.mkdir(parents=True, exist_ok=True)
rng = random.Random(FEM["seed"])
righe_metadata = []
for indice in range(1, FEM["num_run"] + 1):
id_run = f"run_{indice:04d}"
print(f"Integrazione FEM di {id_run} in corso...")
riga = esporta_run(id_run, cartella / f"{id_run}.csv", mesh_dati, rng)
print(
f" T vera max al sensore: {riga['T_vera_max_sensore_C']:.1f} °C — "
f"T misurata max: {riga['T_misurata_max_sensore_C']:.1f} °C"
)
righe_metadata.append(riga)
percorso_metadata = cartella / "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"Esportate {len(righe_metadata)} analisi in: {cartella.resolve()}")
print(f"Metadata: {percorso_metadata.resolve()}")
if __name__ == "__main__":
main()
+460
View File
@@ -0,0 +1,460 @@
# Solutore termico transitorio a elementi finiti sulla mesh shell della
# fascetta (vedi mesh.py per la geometria).
#
# Il dominio è la superficie media del cilindro, con x e la coordinata
# circonferenziale s = R·theta entrambe risolte, mentre lo spessore è
# collassato (temperatura uniforme attraverso la parete).
#
# Il collasso in spessore è lecito: con spessore 0.18 mm il numero di Biot
# h·t/k vale ~1e-7 e il tempo di diffusione attraverso la parete t²/alpha è
# di pochi millisecondi, molto più rapido del transito delle sorgenti. Di
# conseguenza la profondità di penetrazione del riscaldamento a induzione non
# entra nel modello: conta solo il flusso totale assorbito per unità di
# superficie. La skin depth è comunque calcolata dalla frequenza dell'induttore
# come grandezza diagnostica (`calcola_skin_depth_m`).
#
# Equazione risolta (per unità di superficie media):
#
# rho·cp·t ∂T/∂t = k·t ∇²T + q(x, s, t) - (h_est + h_int)·(T - T_amb)
#
# con convezione aggiuntiva sui due bordi anulari x = 0 e x = L (coefficiente
# h_bordi su un'area pari a spessore × perimetro). L'impronta della sorgente
# è una gaussiana isotropa nel piano (x, s), con la distanza circonferenziale
# valutata sull'immagine più vicina perché la superficie è chiusa.
#
# Discretizzazione: elementi shell quadrangolari bilineari a 4 nodi. Gli
# elementi della mesh sono tutti rettangoli identici dx × ds, quindi le
# matrici di elemento sono calcolate una volta sola e assemblate in forma
# vettorizzata. Integrazione temporale con Eulero implicito: la matrice di
# sistema è costante e viene fattorizzata LU una volta per run.
import math
import random
from copy import deepcopy
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import splu
from config import ARIA, FASCETTA, FEM, RANDOMIZZAZIONE, SENSORE, SORGENTE
from materials import MATERIALI
from mesh import genera_mesh
MU0 = 4.0 * math.pi * 1e-7
# Matrici di riferimento dell'elemento rettangolare bilineare a 4 nodi, con
# nodi locali in ordine antiorario: (0,0), (a,0), (a,b), (0,b).
#
# MASSA_RIF va moltiplicata per l'area a·b e dà l'integrale di N_i·N_j;
# RIGIDEZZA_X per b/a e RIGIDEZZA_S per a/b, e la loro somma dà l'integrale
# di grad(N_i)·grad(N_j).
MASSA_RIF = np.array(
[[4.0, 2.0, 1.0, 2.0],
[2.0, 4.0, 2.0, 1.0],
[1.0, 2.0, 4.0, 2.0],
[2.0, 1.0, 2.0, 4.0]]
) / 36.0
RIGIDEZZA_X = np.array(
[[2.0, -2.0, -1.0, 1.0],
[-2.0, 2.0, 1.0, -1.0],
[-1.0, 1.0, 2.0, -2.0],
[1.0, -1.0, -2.0, 2.0]]
) / 6.0
RIGIDEZZA_S = np.array(
[[2.0, 1.0, -1.0, -2.0],
[1.0, 2.0, -2.0, -1.0],
[-1.0, -2.0, 2.0, 1.0],
[-2.0, -1.0, 1.0, 2.0]]
) / 6.0
def calcola_skin_depth_m(materiale: dict, frequenza_hz: float) -> float:
# Skin depth elettromagnetica approssimata:
# delta = sqrt(2 * rho_e / (omega * mu))
#
# Semplificata. Per acciai ferromagnetici il comportamento reale
# è fortemente non lineare con temperatura e campo magnetico.
#
# Nel modello shell è una grandezza diagnostica: serve a verificare che la
# deposizione del calore resti confinata entro uno spessore confrontabile
# con la parete, non entra nell'equazione risolta.
rho_e = materiale["resistivita_elettrica_ohm_m"]
mu_r = materiale["permeabilita_relativa"]
omega = 2.0 * math.pi * frequenza_hz
mu = MU0 * mu_r
return math.sqrt(2.0 * rho_e / (omega * mu))
def quantizza(valore: float, passo: float) -> float:
if passo <= 0.0:
return valore
return round(valore / passo) * passo
def configurazione_randomizzata(rng: random.Random | None = None) -> dict:
"""Copia dei dizionari di configurazione usati da un'analisi FEM.
Se `RANDOMIZZAZIONE["abilitata"]`, velocità, flusso di picco, sigma dello
spot, offset y del percorso, temperatura ambiente e rumore del sensore sono
perturbati con estrazioni dall'RNG passato: analisi diverse dello stesso
setup differiscono come differirebbero due passaggi reali.
"""
fascetta = deepcopy(FASCETTA)
aria = deepcopy(ARIA)
sorgente = deepcopy(SORGENTE)
sensore = deepcopy(SENSORE)
if RANDOMIZZAZIONE.get("abilitata", False):
if rng is None:
rng = random.Random(FEM["seed"])
def perturba_rel(valore: float, std_rel: float, fattore_min: float = 0.1) -> float:
fattore = rng.gauss(1.0, std_rel)
fattore = max(fattore_min, fattore)
return valore * fattore
sorgente["velocita_m_s"] = perturba_rel(
sorgente["velocita_m_s"],
RANDOMIZZAZIONE["velocita_std_rel"],
)
sorgente["flusso_termico_picco_W_m2"] = perturba_rel(
sorgente["flusso_termico_picco_W_m2"],
RANDOMIZZAZIONE["flusso_picco_std_rel"],
)
sorgente["sigma_punto_m"] = perturba_rel(
sorgente["sigma_punto_m"],
RANDOMIZZAZIONE["sigma_punto_std_rel"],
)
sorgente["offset_y_percorso_m"] = rng.uniform(
-RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
)
aria["temperatura_ambiente_C"] += rng.gauss(
0.0,
RANDOMIZZAZIONE["temperatura_ambiente_std_C"],
)
sensore["rumore_std_C"] = perturba_rel(
sensore["rumore_std_C"],
RANDOMIZZAZIONE["rumore_sensore_std_rel"],
fattore_min=0.0,
)
return {
"fascetta": fascetta,
"aria": aria,
"sorgente": sorgente,
"sensore": sensore,
}
def _spread_sorgenti_m(sorgente: dict) -> float:
# Distanza lungo x tra la prima e l'ultima sorgente del gruppo.
numero_sorgenti = sorgente.get("numero_sorgenti", 1)
distanza = sorgente.get("distanza_sorgenti_m", 0.0)
return (numero_sorgenti - 1) * distanza
def _x_riferimento_iniziale_m(sorgente: dict, x_sensore_m: float) -> float:
# Posizione a t=0 della sorgente di indice 0. Le sorgenti i sono a
# x_i = riferimento + i * distanza, quindi per v >= 0 l'indice 0 è la
# più arretrata nel verso di marcia, per v < 0 è la più avanzata.
# x_inizio_m è la distanza dal sensore della sorgente più avanzata
# (quella che lo raggiunge per prima).
spread = _spread_sorgenti_m(sorgente)
x_inizio = sorgente["x_inizio_m"]
if sorgente["velocita_m_s"] >= 0:
return (x_sensore_m - x_inizio) - spread
return x_sensore_m + x_inizio
def _x_riferimento_finale_m(sorgente: dict, x_sensore_m: float) -> float:
# Posizione di fine corsa della sorgente di indice 0. x_fine_m è la
# distanza dal sensore della sorgente più arretrata nel verso di marcia
# (quella che lo supera per ultima).
spread = _spread_sorgenti_m(sorgente)
x_fine = sorgente["x_fine_m"]
if sorgente["velocita_m_s"] >= 0:
return x_sensore_m + x_fine
return (x_sensore_m - x_fine) - spread
def _intervallo_attivo(inizio: float, fine: float, v: float, x_m: float) -> bool:
if v >= 0:
return inizio <= x_m <= fine
return fine <= x_m <= inizio
def _assembla(elementi: np.ndarray, matrice_elemento: np.ndarray, n_nodi: int):
# Assembla una matrice globale sparsa a partire da un'unica matrice di
# elemento 4×4, uguale per tutti gli elementi. I contributi ripetuti sullo
# stesso (riga, colonna) vengono sommati dal formato COO.
n_elementi = len(elementi)
righe = np.broadcast_to(elementi[:, :, None], (n_elementi, 4, 4)).ravel()
colonne = np.broadcast_to(elementi[:, None, :], (n_elementi, 4, 4)).ravel()
valori = np.broadcast_to(matrice_elemento, (n_elementi, 4, 4)).ravel()
return sp.coo_matrix(
(valori, (righe, colonne)), shape=(n_nodi, n_nodi)
).tocsr()
def _massa_anello(indici: np.ndarray, lunghezza_segmento_m: float, n_nodi: int):
# Matrice di massa 1D su un anello chiuso di nodi (integrale di N_i·N_j
# lungo la circonferenza), usata per la convezione sui bordi x = 0 e x = L.
n = len(indici)
successivi = indici[(np.arange(n) + 1) % n]
coefficiente = lunghezza_segmento_m / 6.0
righe = np.concatenate([indici, indici, successivi, successivi])
colonne = np.concatenate([indici, successivi, indici, successivi])
valori = coefficiente * np.concatenate(
[np.full(n, 2.0), np.full(n, 1.0), np.full(n, 1.0), np.full(n, 2.0)]
)
return sp.coo_matrix(
(valori, (righe, colonne)), shape=(n_nodi, n_nodi)
).tocsr()
def prepara_stato_fem(
fascetta: dict,
aria: dict,
sorgente: dict,
sensore: dict,
mesh_dati: dict | None = None,
dt_s: float | None = None,
) -> dict:
"""Assembla le matrici FEM e fattorizza il sistema implicito.
Restituisce lo stato costante per l'integrazione temporale: matrici,
solutore LU, coordinate nodali, peso circonferenziale della sorgente e
indice del nodo osservato dal sensore.
"""
if mesh_dati is None:
mesh_dati = genera_mesh(fascetta)
if dt_s is None:
dt_s = FEM["dt_s"]
materiale = MATERIALI[fascetta["materiale"]]
if sorgente["skin_depth_fissa_m"] is None:
skin_depth_m = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
else:
skin_depth_m = float(sorgente["skin_depth_fissa_m"])
k = materiale["conducibilita_termica_W_mK"]
rho = materiale["densita_kg_m3"]
cp = materiale["calore_specifico_J_kgK"]
spessore_m = mesh_dati["spessore_m"]
raggio_m = mesh_dati["raggio_m"]
circonferenza_m = 2.0 * math.pi * raggio_m
n_x = mesh_dati["n_elementi_x"]
n_theta = mesh_dati["n_elementi_circonferenza"]
elementi = mesh_dati["elementi"]
indice_nodo = mesh_dati["indice_nodo"]
n_nodi = len(mesh_dati["nodi"])
dx_m = mesh_dati["lunghezza_m"] / n_x
ds_m = circonferenza_m / n_theta
h_esterno = aria["h_esterno_W_m2K"]
h_interno = aria["h_interno_W_m2K"]
h_bordi = aria["h_bordi_W_m2K"]
T_ambiente = aria["temperatura_ambiente_C"]
# Matrice di massa di superficie: riusata per la capacità termica, per la
# convezione sulle due facce e per il carico della sorgente.
M_area = _assembla(elementi, MASSA_RIF * (dx_m * ds_m), n_nodi)
K_cond = _assembla(
elementi,
k * spessore_m * (RIGIDEZZA_X * (ds_m / dx_m) + RIGIDEZZA_S * (dx_m / ds_m)),
n_nodi,
)
# Convezione sui bordi anulari: l'area di scambio di ogni segmento è
# spessore × lunghezza del segmento circonferenziale.
H_bordi = h_bordi * spessore_m * (
_massa_anello(indice_nodo[0, :], ds_m, n_nodi)
+ _massa_anello(indice_nodo[-1, :], ds_m, n_nodi)
)
C_su_dt = (rho * cp * spessore_m / dt_s) * M_area
A = C_su_dt + K_cond + (h_esterno + h_interno) * M_area + H_bordi
# Termine noto costante della convezione: (somma delle righe) × h × T_amb,
# perché la somma delle righe della matrice di massa è l'integrale di N_i.
area_nodale = np.asarray(M_area.sum(axis=1)).ravel()
bordo_nodale = np.asarray(H_bordi.sum(axis=1)).ravel()
carico_convezione = (h_esterno + h_interno) * area_nodale * T_ambiente
carico_convezione += bordo_nodale * T_ambiente
# Peso circonferenziale dell'impronta della sorgente: dipende solo
# dall'offset y del percorso e dalla sigma dello spot, entrambi costanti
# durante il run, quindi è valutato una volta sola qui. La distanza è
# quella dell'immagine più vicina, perché la superficie è chiusa.
sigma_m = sorgente["sigma_punto_m"]
ds_sorgente = raggio_m * mesh_dati["theta_nodi_rad"] - sorgente["offset_y_percorso_m"]
ds_sorgente -= circonferenza_m * np.round(ds_sorgente / circonferenza_m)
peso_circonferenziale = np.exp(-0.5 * (ds_sorgente ** 2) / (sigma_m * sigma_m))
# Nodo osservato dal sensore: x più vicino a x_mm sul piano theta = 0.
x_sensore_m = sensore["x_mm"] / 1000.0
i_sensore = int(np.argmin(np.abs(mesh_dati["x_nodi_m"] - x_sensore_m)))
indice_sensore = int(indice_nodo[i_sensore, 0])
return {
"mesh": mesh_dati,
"dt_s": dt_s,
"M_area": M_area,
"C_su_dt": C_su_dt,
"carico_convezione": carico_convezione,
"solutore": splu(sp.csc_matrix(A)),
"x_nodi_m": mesh_dati["x_nodi_m"],
"peso_circonferenziale": peso_circonferenziale,
"forma_griglia": (n_x + 1, n_theta),
"circonferenza_m": circonferenza_m,
"area_nodale_m2": area_nodale,
"T_ambiente_C": T_ambiente,
"x_sensore_m": x_sensore_m,
"indice_sensore": indice_sensore,
"n_nodi": n_nodi,
# Diagnostica: non entra nell'equazione risolta.
"skin_depth_m": skin_depth_m,
}
def flusso_nodale_W_m2(
sorgente: dict,
stato: dict,
t_s: float,
) -> tuple[float, np.ndarray]:
"""Flusso termico assorbito nei nodi all'istante t_s.
Restituisce (x_riferimento_m, q_nodale) dove x_riferimento_m è la
posizione della sorgente di indice 0. Ogni sorgente attiva contribuisce
con un'impronta gaussiana isotropa in (x, s).
L'impronta gaussiana isotropa è separabile, q(x, s) = q_x(x) · q_s(s), e
il fattore circonferenziale q_s è costante nel tempo (precalcolato in
`prepara_stato_fem`). Qui si valuta quindi solo il fattore assiale sui
n_x + 1 nodi distinti in x e si espande con un prodotto esterno, invece di
valutare l'esponenziale su tutti i nodi della mesh.
"""
x_sensore_m = stato["x_sensore_m"]
x_rif_iniziale = _x_riferimento_iniziale_m(sorgente, x_sensore_m)
x_rif_finale = _x_riferimento_finale_m(sorgente, x_sensore_m)
x_riferimento = x_rif_iniziale + sorgente["velocita_m_s"] * t_s
numero_sorgenti = sorgente.get("numero_sorgenti", 1)
distanza = sorgente.get("distanza_sorgenti_m", 0.0)
v = sorgente["velocita_m_s"]
zero_dopo_fine = sorgente.get("zero_dopo_fine", True)
sigma = sorgente["sigma_punto_m"]
q_picco = sorgente["flusso_termico_picco_W_m2"] * sorgente["efficienza_riscaldamento"]
peso_assiale = np.zeros(len(stato["x_nodi_m"]))
for i in range(numero_sorgenti):
x_i = x_riferimento + i * distanza
if zero_dopo_fine:
inizio_i = x_rif_iniziale + i * distanza
fine_i = x_rif_finale + i * distanza
if not _intervallo_attivo(inizio_i, fine_i, v, x_i):
continue
dx = stato["x_nodi_m"] - x_i
peso_assiale += np.exp(-0.5 * (dx * dx) / (sigma * sigma))
q_nodale = q_picco * np.outer(peso_assiale, stato["peso_circonferenziale"])
return x_riferimento, q_nodale.ravel()
def passo_implicito_fem(stato: dict, T: np.ndarray, q_nodale: np.ndarray) -> np.ndarray:
# Avanza il campo nodale di un passo dt: termine noto = capacità sul campo
# precedente + carico della sorgente + carico di convezione, poi risoluzione
# del sistema già fattorizzato.
rhs = stato["C_su_dt"] @ T
rhs += stato["M_area"] @ q_nodale
rhs += stato["carico_convezione"]
return stato["solutore"].solve(rhs)
def campo_iniziale(stato: dict) -> np.ndarray:
return np.full(stato["n_nodi"], stato["T_ambiente_C"], dtype=float)
def simula_campo_fem(
dt_frame_s: float,
cfg: dict | None = None,
durata_s: float | None = None,
mesh_dati: dict | None = None,
rng: random.Random | None = None,
salva_campi: bool = True,
) -> dict:
"""Integra il campo FEM fino a durata_s, campionando 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.
Con `salva_campi = False` i campi nodali non vengono accumulati (utile per
l'export CSV, che usa solo le serie scalari): `campi` resta una lista vuota.
"""
if rng is None:
rng = random.Random(FEM["seed"])
if cfg is None:
cfg = configurazione_randomizzata(rng)
if durata_s is None:
durata_s = FEM["durata_s"]
fascetta = cfg["fascetta"]
aria = cfg["aria"]
sorgente = cfg["sorgente"]
sensore = cfg["sensore"]
stato = prepara_stato_fem(fascetta, aria, sorgente, sensore, mesh_dati=mesh_dati)
dt = stato["dt_s"]
indice_sensore = stato["indice_sensore"]
T = campo_iniziale(stato)
T_sensore = float(T[indice_sensore])
tau_sensore = max(sensore["costante_tempo_s"], 1e-9)
tempi, campi, x_riferimenti = [], [], []
T_vere, T_lette, T_massime = [], [], []
prossimo_frame_t = 0.0
t = 0.0
while t <= durata_s + 1e-12:
x_rif, q_nodale = flusso_nodale_W_m2(sorgente, stato, t)
T = passo_implicito_fem(stato, T, q_nodale)
T_sensore += (T[indice_sensore] - T_sensore) * dt / tau_sensore
if t + 1e-12 >= prossimo_frame_t:
letta = T_sensore + rng.gauss(0.0, sensore["rumore_std_C"])
letta = quantizza(letta, sensore["quantizzazione_C"])
tempi.append(t)
if salva_campi:
campi.append(T.copy())
x_riferimenti.append(x_rif)
T_vere.append(float(T[indice_sensore]))
T_lette.append(letta)
T_massime.append(float(T.max()))
prossimo_frame_t += dt_frame_s
t += dt
return {
"stato": stato,
"mesh": stato["mesh"],
"tempi": np.array(tempi),
"campi": campi,
"x_riferimenti": np.array(x_riferimenti),
"T_vere": np.array(T_vere),
"T_lette": np.array(T_lette),
"T_massime": np.array(T_massime),
"cfg": cfg,
"sorgente": sorgente,
"T_ambiente_C": stato["T_ambiente_C"],
"skin_depth_m": stato["skin_depth_m"],
}
+61
View File
@@ -0,0 +1,61 @@
# Punto di ingresso unico del progetto: seleziona da riga di comando quale
# azione eseguire.
#
# python main.py mesh # disegna la sola mesh a elementi shell
# python main.py fem # animazione 3D del campo FEM sulla mesh shell
# python main.py csv # esporta in CSV la serie temporale del FEM
#
# Senza argomenti stampa l'elenco delle azioni disponibili.
import argparse
import sys
def _azione_mesh() -> None:
import plot_mesh
plot_mesh.main()
def _azione_fem() -> None:
import plot_animazione_fem
plot_animazione_fem.main()
def _azione_csv() -> None:
import esporta_csv
esporta_csv.main()
# Chiave da riga di comando -> (funzione, descrizione mostrata nell'help).
AZIONI = {
"mesh": (_azione_mesh, "Disegna la mesh a elementi shell della fascetta"),
"fem": (_azione_fem, "Animazione 3D del campo FEM sulla mesh a elementi shell"),
"csv": (_azione_csv, "Esporta in CSV la temperatura simulata e la lettura del sensore"),
}
def main() -> None:
parser = argparse.ArgumentParser(
description="Analisi termica FEM della fascetta riscaldata a induzione.",
formatter_class=argparse.RawTextHelpFormatter,
)
parser.add_argument(
"azione",
nargs="?",
choices=list(AZIONI),
help="\n".join(f"{nome}: {descrizione}" for nome, (_, descrizione) in AZIONI.items()),
)
argomenti = parser.parse_args()
if argomenti.azione is None:
parser.print_help()
sys.exit(1)
AZIONI[argomenti.azione][0]()
if __name__ == "__main__":
main()
+1 -1
View File
@@ -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
+114
View File
@@ -0,0 +1,114 @@
# Generazione della mesh a elementi shell quadrangolari della fascetta.
#
# La fascetta è modellata come la superficie media di un cilindro: gli
# elementi sono quadrilateri a 4 nodi disposti su una griglia strutturata
# n_elementi_x × n_elementi_circonferenza. Lo spessore della parete non è
# discretizzato, è un attributo degli elementi shell.
#
# Sistema di coordinate globale (coerente con plot_animazione_3d.py):
# - x = asse del cilindro, da 0 a lunghezza
# - y, z = piano della sezione circolare, con y = R·sin(theta), z = R·cos(theta)
# - theta = 0 corrisponde al piano del sensore (y = 0), cresce in verso
# antiorario nel piano y-z
#
# La mesh è chiusa lungo la circonferenza: l'ultima colonna di elementi
# richiude sulla prima riga di nodi, senza nodi duplicati sulla cucitura.
import numpy as np
from config import FASCETTA, MESH
def raggio_medio_m(fascetta: dict = FASCETTA) -> float:
# Raggio della superficie media della parete: il diametro configurato è
# quello esterno, la shell sta a metà dello spessore.
diametro_m = fascetta["diametro_mm"] / 1000.0
spessore_m = fascetta["spessore_mm"] / 1000.0
return (diametro_m - spessore_m) / 2.0
def genera_mesh(fascetta: dict = FASCETTA, mesh: dict = MESH) -> dict:
"""Costruisce la mesh shell quadrangolare della fascetta.
Restituisce un dizionario con:
- nodi: array (n_nodi, 3) con le coordinate globali [m]
- elementi: array (n_elementi, 4) con gli indici dei nodi, in verso
antiorario visto dall'esterno del cilindro
- indice_nodo: array (n_x + 1, n_theta) che mappa (i, j) sull'indice
globale del nodo
- x_nodi_m, theta_nodi_rad: coordinate parametriche della griglia
- spessore_m, raggio_m, lunghezza_m: geometria della shell
"""
n_x = int(mesh["n_elementi_x"])
n_theta = int(mesh["n_elementi_circonferenza"])
if n_x < 1 or n_theta < 3:
raise ValueError(
"Servono almeno 1 elemento lungo x e 3 lungo la circonferenza."
)
lunghezza_m = fascetta["lunghezza_mm"] / 1000.0
spessore_m = fascetta["spessore_mm"] / 1000.0
raggio_m = raggio_medio_m(fascetta)
# Lungo x la griglia è aperta (n_x + 1 file di nodi), lungo theta è
# chiusa (n_theta file, l'ultima si ricongiunge alla prima).
x_nodi_m = np.linspace(0.0, lunghezza_m, n_x + 1)
theta_nodi_rad = np.linspace(0.0, 2.0 * np.pi, n_theta, endpoint=False)
X, Theta = np.meshgrid(x_nodi_m, theta_nodi_rad, indexing="ij")
nodi = np.column_stack(
(
X.ravel(),
(raggio_m * np.sin(Theta)).ravel(),
(raggio_m * np.cos(Theta)).ravel(),
)
)
indice_nodo = np.arange((n_x + 1) * n_theta).reshape(n_x + 1, n_theta)
i = np.arange(n_x)[:, None]
j = np.arange(n_theta)[None, :]
j_succ = (j + 1) % n_theta
elementi = np.stack(
(
np.broadcast_to(indice_nodo[i, j], (n_x, n_theta)),
np.broadcast_to(indice_nodo[i + 1, j], (n_x, n_theta)),
np.broadcast_to(indice_nodo[i + 1, j_succ], (n_x, n_theta)),
np.broadcast_to(indice_nodo[i, j_succ], (n_x, n_theta)),
),
axis=-1,
).reshape(-1, 4)
return {
"nodi": nodi,
"elementi": elementi,
"indice_nodo": indice_nodo,
"x_nodi_m": x_nodi_m,
"theta_nodi_rad": theta_nodi_rad,
"n_elementi_x": n_x,
"n_elementi_circonferenza": n_theta,
"spessore_m": spessore_m,
"raggio_m": raggio_m,
"lunghezza_m": lunghezza_m,
}
def riepilogo_mesh(mesh_dati: dict) -> str:
n_nodi = len(mesh_dati["nodi"])
n_elementi = len(mesh_dati["elementi"])
passo_x_mm = 1000.0 * mesh_dati["lunghezza_m"] / mesh_dati["n_elementi_x"]
passo_circ_mm = (
1000.0
* 2.0
* np.pi
* mesh_dati["raggio_m"]
/ mesh_dati["n_elementi_circonferenza"]
)
return (
f"Mesh shell: {n_elementi} elementi quadrangolari, {n_nodi} nodi\n"
f" griglia: {mesh_dati['n_elementi_x']} (x) × "
f"{mesh_dati['n_elementi_circonferenza']} (circonferenza)\n"
f" passo elemento: {passo_x_mm:.2f} mm (x), {passo_circ_mm:.2f} mm (circonferenza)\n"
f" raggio medio: {1000.0 * mesh_dati['raggio_m']:.2f} mm, "
f"spessore shell: {1000.0 * mesh_dati['spessore_m']:.3f} mm"
)
-230
View File
@@ -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()
+198
View File
@@ -0,0 +1,198 @@
# Animazione 3D isometrica della fascetta con il campo di temperatura calcolato
# a elementi finiti sulla mesh shell (fem.py).
#
# 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.
from pathlib import Path
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
from matplotlib.animation import FuncAnimation, PillowWriter
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from config import USCITA
from fem import simula_campo_fem
from mesh import genera_mesh, riepilogo_mesh
# 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
# Millisecondi tra i fotogrammi in riproduzione.
INTERVALLO_RIPRODUZIONE_MS = 60
# Direzione della luce (x = asse del cilindro, y-z = piano della sezione) e
# quota di luce ambiente, per dare volume alla superficie.
DIREZIONE_LUCE = np.array([-0.25, 0.45, 0.85])
LUCE_AMBIENTE = 0.45
# Se True disegna il reticolo degli elementi sopra il campo. Con la mesh fitta
# i bordi degli elementi sul lato nascosto rendono la superficie confusa,
# quindi per default i bordi prendono il colore della faccia: così le facce
# adiacenti combaciano senza lasciare fessure di antialiasing.
MOSTRA_BORDI_ELEMENTI = False
def _illuminazione(mesh_dati: dict) -> np.ndarray:
# Fattore moltiplicativo di luminosità per elemento, dalla normale
# radiale uscente valutata nel baricentro dell'elemento.
facce = mesh_dati["nodi"][mesh_dati["elementi"]]
baricentri = facce.mean(axis=1)
normali = baricentri.copy()
# La normale della superficie cilindrica è radiale: nessuna componente
# lungo l'asse x.
normali[:, 0] = 0.0
normali /= np.linalg.norm(normali, axis=1, keepdims=True)
luce = DIREZIONE_LUCE / np.linalg.norm(DIREZIONE_LUCE)
diffusa = np.clip(normali @ luce, 0.0, None)
return LUCE_AMBIENTE + (1.0 - LUCE_AMBIENTE) * diffusa
def main() -> None:
mesh_dati = genera_mesh()
print(riepilogo_mesh(mesh_dati))
print("Integrazione FEM in corso...")
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))
indici_frame = list(range(indice_inizio, len(tempi)))
T_ambiente = dati["T_ambiente_C"]
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\n"
f"Skin depth (diagnostica): {1000.0 * dati['skin_depth_m']:.3f} mm"
)
elementi = mesh_dati["elementi"]
facce = mesh_dati["nodi"][elementi]
illuminazione = _illuminazione(mesh_dati)[:, None]
raggio_m = mesh_dati["raggio_m"]
lunghezza_m = mesh_dati["lunghezza_m"]
sorgente = dati["sorgente"]
numero_sorgenti = sorgente.get("numero_sorgenti", 1)
distanza_m = sorgente.get("distanza_sorgenti_m", 0.0)
theta_sorgente = sorgente["offset_y_percorso_m"] / raggio_m
norm = matplotlib.colors.Normalize(vmin=T_ambiente, vmax=T_max)
cmap = matplotlib.colormaps["inferno"]
fig = plt.figure(figsize=(9, 8.5))
griglia = fig.add_gridspec(2, 1, height_ratios=[2.6, 1.0], hspace=0.05)
ax = fig.add_subplot(griglia[0], projection="3d")
ax_storia = fig.add_subplot(griglia[1])
ax.view_init(elev=35.264, azim=45)
# zoom > 1 riempie il riquadro: con set_axis_off gli assi 3D lascerebbero
# altrimenti molto margine vuoto sopra il pannello della storia.
ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m), zoom=1.12)
ax.set_axis_off()
ax.set_xlim(0.0, lunghezza_m)
ax.set_ylim(-raggio_m, raggio_m)
ax.set_zlim(-raggio_m, raggio_m)
collezione = Poly3DCollection(facce, linewidths=0.25, shade=False)
ax.add_collection3d(collezione)
marker_sorgenti = ax.plot(
[], [], [], "o", color="cyan", markersize=6, zorder=10
)[0]
mappabile = cm.ScalarMappable(cmap=cmap, norm=norm)
mappabile.set_array([])
# La colorbar è agganciata a entrambi i pannelli per non restringere solo
# quello 3D: così i due riquadri restano della stessa larghezza.
barra = fig.colorbar(mappabile, ax=(ax, ax_storia), shrink=0.6, pad=0.02, aspect=30)
barra.set_label("T [°C]")
titolo = ax.set_title("")
# Pannello inferiore: storia della temperatura nel punto osservato dal
# sensore. T_vere è il valore nodale calcolato dal FEM, T_lette lo stesso
# 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.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])
ax_storia.set_ylim(
T_ambiente - 0.05 * (T_max - T_ambiente),
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)
def disegna_frame(k: int):
T_elementi = dati["campi"][k][elementi].mean(axis=1)
colori = cmap(norm(T_elementi))
# Ombreggiatura moltiplicativa sui soli canali RGB.
colori[:, :3] *= illuminazione
collezione.set_facecolor(colori)
collezione.set_edgecolor("#40404060" if MOSTRA_BORDI_ELEMENTI else colori)
x_sorgenti_m = dati["x_riferimenti"][k] + np.arange(numero_sorgenti) * distanza_m
visibili = (x_sorgenti_m >= 0.0) & (x_sorgenti_m <= lunghezza_m)
n_visibili = int(visibili.sum())
marker_sorgenti.set_data_3d(
x_sorgenti_m[visibili],
np.full(n_visibili, raggio_m * np.sin(theta_sorgente) * 1.06),
np.full(n_visibili, raggio_m * np.cos(theta_sorgente) * 1.06),
)
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([tempi[k], tempi[k]])
titolo.set_text(
f"FEM shell — t = {tempi[k]:.2f} s — "
f"T sensore = {dati['T_lette'][k]:.1f} °C"
)
return collezione, marker_sorgenti, linea_vera, linea_letta, cursore, titolo
animazione = FuncAnimation(
fig, disegna_frame, frames=indici_frame,
interval=INTERVALLO_RIPRODUZIONE_MS, blit=False,
)
if matplotlib.get_backend().lower() == "agg":
cartella = Path(USCITA["cartella"])
cartella.mkdir(parents=True, exist_ok=True)
percorso = cartella / "animazione_fem.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()
-48
View File
@@ -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()
+60
View File
@@ -0,0 +1,60 @@
# Visualizzazione della sola mesh a elementi shell della fascetta, in vista
# isometrica: facce chiare con i bordi degli elementi in evidenza.
from pathlib import Path
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
def main() -> None:
dati = genera_mesh()
print(riepilogo_mesh(dati))
nodi = dati["nodi"]
facce = nodi[dati["elementi"]]
raggio_m = dati["raggio_m"]
lunghezza_m = dati["lunghezza_m"]
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_axis_off()
collezione = Poly3DCollection(
facce,
facecolors="#d9e3f0",
edgecolors="#2b3a55",
linewidths=0.4,
shade=True,
lightsource=matplotlib.colors.LightSource(azdeg=315, altdeg=45),
)
ax.add_collection3d(collezione)
ax.set_xlim(0.0, lunghezza_m)
ax.set_ylim(-raggio_m, raggio_m)
ax.set_zlim(-raggio_m, raggio_m)
ax.set_title(
f"Mesh shell: {dati['n_elementi_x']} × "
f"{dati['n_elementi_circonferenza']} elementi quadrangolari"
)
if matplotlib.get_backend().lower() == "agg":
cartella = Path(USCITA["cartella"])
cartella.mkdir(parents=True, exist_ok=True)
percorso = cartella / "mesh.png"
fig.savefig(percorso, dpi=150)
print(f"Backend non interattivo: immagine salvata in {percorso}")
return
plt.show()
if __name__ == "__main__":
main()
-470
View File
@@ -1,470 +0,0 @@
import csv
import math
import random
import shutil
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 main() -> None:
cartella_output = Path(SIMULAZIONE["cartella_output"])
if cartella_output.exists():
shutil.rmtree(cartella_output)
cartella_output.mkdir(parents=True, exist_ok=True)
rng = random.Random(SIMULAZIONE["seed"])
righe_metadata = []
for i in range(1, SIMULAZIONE["num_run"] + 1):
cfg_run = configurazione_randomizzata(i, rng)
percorso_csv = cartella_output / f"{cfg_run['id_run']}.csv"
righe_metadata.append(simula_singolo(cfg_run, percorso_csv, rng))
percorso_metadata = cartella_output / "metadata.csv"
with percorso_metadata.open("w", newline="") as f:
writer = csv.DictWriter(f, fieldnames=list(righe_metadata[0].keys()))
writer.writeheader()
writer.writerows(righe_metadata)
print(f"Generati {len(righe_metadata)} run in: {cartella_output.resolve()}")
print(f"Metadata: {percorso_metadata.resolve()}")
if __name__ == "__main__":
main()