From 7883f3deadb9f88dd81463ef9bcd719d09583176 Mon Sep 17 00:00:00 2001 From: Davide Grilli Date: Mon, 3 Aug 2026 00:18:00 +0200 Subject: [PATCH] temp --- CLAUDE.md | 50 ++++-- README.md | 293 ++++++++++++++++++++++--------- assemblaggio.py | 151 ++++++++++++++++ config.py | 79 +++++---- elementi_shell.py | 135 +++++++++++++++ mesh.py | 119 +++++++++++++ plot_animazione.py | 169 +++++++++++------- plot_animazione_3d.py | 84 ++++----- sensore.py | 60 +++++++ simulate.py | 392 ++++++++++++++++-------------------------- solutore.py | 39 +++++ sorgente.py | 125 ++++++++++++++ valida_shell.py | 252 +++++++++++++++++++++++++++ 13 files changed, 1460 insertions(+), 488 deletions(-) create mode 100644 assemblaggio.py create mode 100644 elementi_shell.py create mode 100644 mesh.py create mode 100644 sensore.py create mode 100644 solutore.py create mode 100644 sorgente.py create mode 100644 valida_shell.py diff --git a/CLAUDE.md b/CLAUDE.md index 8770517..e29c4bf 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -23,40 +23,59 @@ python simulate.py # Visualizzare il primo run python plot_csv.py + +# Animazione della superficie sviluppata e vista 3D isometrica +python plot_animazione.py +python plot_animazione_3d.py ``` Attivare sempre il venv (`source .venv/bin/activate`) prima di eseguire qualsiasi comando Python. -Non sono configurati test o linter. +Non sono configurati linter. La validazione del solutore si esegue a mano: + +```bash +python valida_shell.py +``` + +Va rieseguita dopo ogni modifica alla fisica o all'assemblaggio: confronta le matrici e il solutore con soluzioni analitiche (raffreddamento convettivo, modi sinusoidali circonferenziali, conservazione dell'energia, potenza della sorgente, periodicità) e stampa gli errori. ## Architettura Generatore di dataset per misurazioni termiche pseudo-realistiche 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 medio, spessore e lunghezza configurabili. Il dominio simulato è l'intera superficie cilindrica, parametrizzata da x (lunghezza, le sorgenti si muovono in direzione -x) e theta (circonferenza, periodica). Lo spessore non è discretizzato: la parete è termicamente sottile, quindi faccia esterna e faccia interna hanno la stessa temperatura. Il sensore è un pirometro a infrarossi dentro la fascetta che misura la superficie interna in un punto x fisso, sulla generatrice theta = 0; l'offset circonferenziale delle sorgenti è misurato rispetto a quella generatrice. **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 +3. `mesh.py` — nodi, connettività Q4 e periodicità circonferenziale della shell cilindrica +4. `elementi_shell.py` — funzioni di forma bilineari, quadratura di Gauss 2×2, matrici del singolo elemento +5. `assemblaggio.py` — assemblaggio delle matrici globali e del vettore della sorgente +6. `sorgente.py` — moto del gruppo di sorgenti e distribuzione gaussiana del flusso +7. `solutore.py` — Eulero implicito e fattorizzazione LU sparsa +8. `sensore.py` — interpolazione della temperatura, inerzia, rumore, quantizzazione +9. `simulate.py` — motore principale: genera N run randomizzati, scrive i CSV, scrive `metadata.csv` +10. `plot_csv.py`, `plot_animazione.py`, `plot_animazione_3d.py` — visualizzazioni autonome -**Pipeline fisica dentro `simula_singolo()` in [simulate.py](simulate.py):** +**Pipeline fisica in [simulate.py](simulate.py):** -- La skin depth è calcolata dalla resistività elettrica del materiale e dalla frequenza di induzione (`calcola_skin_depth_m`) -- Le sorgenti gaussiane in movimento producono un profilo di flusso termico superficiale q(x) sul lato esterno, variabile nel tempo (`profilo_flusso_incidente_W_m2`) -- Quel flusso è ridistribuito volumetricamente attraverso lo spessore con decadimento esponenziale in z (`profilo_deposizione_z_1_m`): q_vol(x, z) = q(x) · p(z) -- Uno schema 2D a volumi finiti con Eulero implicito integra l'equazione del calore su `n_nodi_x × n_nodi_z` celle: `prepara_stato_termico` costruisce griglia, coefficienti e matrice sparsa fattorizzata LU una volta per run (`costruisci_solutore_implicito_2d`, che restituisce l'oggetto `splu`), poi `passo_implicito` avanza il campo risolvendo solo il sistema triangolare -- Le condizioni al contorno sono incorporate nella matrice: convezione su tutti e quattro i lati della sezione, più un termine di conduzione circonferenziale (y) verso il resto della fascetta assunto a temperatura ambiente. Il termine è un'equazione di aletta ricavata sull'intero volume del cilindro: il calore conduce lungo y attraverso l'intero spessore mentre le superfici esterna e interna dell'intero cilindro scambiano per convezione, dando q_y = -(h_esterno + h_interno)/spessore · (T - T_amb), senza parametri di conduzione y configurabili a parte +- `prepara_stato_termico` costruisce la mesh, assembla C, K_cond e K_conv, fattorizza LU una volta per run e prepara l'assemblatore della sorgente e l'interpolatore del sensore +- La conduzione assiale e quella circonferenziale nascono dall'assemblaggio: elementi adiacenti condividono i nodi del bordo comune, quindi non serve alcun accoppiamento esplicito tra elementi +- Convezione su entrambe le facce (h_esterno + h_interno su tutta l'area) e sui due bordi assiali (h_bordi · spessore lungo il perimetro), entrambe incorporate in K_conv e nel termine noto costante `f_ambiente` +- Le sorgenti gaussiane in movimento producono q''(x, theta, t); la distanza circonferenziale è riavvolta su 2πR, quindi la sorgente scalda anche i punti raggiunti dall'altro verso. Il vettore di carico è ∫ N^T q'' dA, integrato con la stessa quadratura di Gauss delle matrici +- Il flusso è separabile in un fattore assiale variabile nel tempo e uno circonferenziale costante: `assembla_sorgente` sfrutta la separazione, non è un'approssimazione +- `passo_termico` avanza il campo risolvendo (C/dt + K) T = C/dt · T_prec + f_ambiente + f_sorgente - La temperatura iniziale del campo è la temperatura ambiente (randomizzata per run) - L'output del sensore aggiunge inerzia del primo ordine, rumore gaussiano e quantizzazione -`prepara_stato_termico` e `passo_implicito` sono condivisi con `plot_animazione.py`, che riproduce la fisica di run_0001 per animare la sezione: ogni modifica alla fisica va fatta lì, non duplicata. +La skin depth (`calcola_skin_depth_m`) e il numero di Fourier sullo spessore (`numero_fourier_spessore`) sono calcolati e registrati in `metadata.csv` come diagnostiche: non entrano nella soluzione, servono a verificare che l'ipotesi di parete termicamente sottile regga. -**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à. +`prepara_stato_termico` e `passo_termico` sono condivisi con `plot_animazione.py`, che riproduce la fisica di run_0001 per animare la superficie, e da lì con `plot_animazione_3d.py`: ogni modifica alla fisica va fatta lì, non duplicata. -**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. +**Randomizzazione per run** (`configurazione_randomizzata`): ogni run perturba velocità, flusso di picco, sigma del punto, offset circonferenziale, 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_max_fascetta_C, T_media_fascetta_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, le dimensioni della mesh, le diagnostiche e le temperature di picco. ## Convenzioni su `config.py` @@ -64,8 +83,11 @@ 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 con una temperatura per nodo: lo spessore non è discretizzato e la temperatura della faccia esterna coincide con quella della faccia interna. L'ipotesi regge finché `numero_fourier_spessore` (α · σ/|v| / spessore²) resta molto maggiore di 1; sotto ~1 servirebbe una formulazione multistrato. Con i valori di default vale circa 5.4, e il gradiente misurato nello spessore col vecchio modello 2D era di 0.002 °C al picco. +- La mesh cilindrica è strutturata e uniforme, quindi tutti gli elementi sono lo stesso rettangolo: le matrici elementari si calcolano una volta e si replicano. Se si introducessero mesh non uniformi o spessore variabile, `assemblaggio.py` dovrebbe ciclare sugli elementi. +- La periodicità circonferenziale è strutturale, non un vincolo imposto: non esiste un nodo separato per theta = 2π, la connettività riavvolge l'ultimo elemento sul primo. - 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. +- La matrice è simmetrica definita positiva: `splu` usa `permc_spec="MMD_AT_PLUS_A"`, che dimezza circa il riempimento rispetto al COLAMD di default. - `simulate.py` cancella e ricrea l'intera cartella di output ad ogni esecuzione (`shutil.rmtree`). - Aggiungere un nuovo materiale richiede solo una nuova voce nel dizionario `MATERIALI` in [materials.py](materials.py); la chiave del materiale va poi impostata in `FASCETTA["materiale"]` in [config.py](config.py). diff --git a/README.md b/README.md index a0c7e71..ced518d 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# Simulatore Termico 2D — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso +# Simulatore Termico a Elementi Shell — 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 @@ -9,105 +9,191 @@ per l'addestramento e la validazione di modelli di stima/regressione termica. 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 medio** (default 70 mm) — il diametro della superficie media, su cui + giace la mesh; +- **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 simulato è l'**intera superficie cilindrica**, discretizzata con elementi +shell quadrilateri. La superficie è parametrizzata da `(x, θ)`: ```text - sorgenti (induttori), in moto verso -x - ▼ ▼ ▼ - (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) + sorgenti (induttori), in moto verso -x + ▼ ▼ ▼ + ╭────────────────────────────────╮ + θ = 0 ─ ─ ─ ┤ ● ← punto osservato dal ├ ─ ─ ─ generatrice del sensore + │ sensore IR (x = 50 mm) │ + ╰────────────────────────────────╯ + x = 0 x = lunghezza + + superficie sviluppata in piano; il bordo superiore e quello inferiore + sono lo stesso punto della fascetta (θ = 0 ≡ θ = 2π) ``` - **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. + direzione `-x` sulla superficie esterna. +- **θ** = coordinata circonferenziale, **periodica** e risolta spazialmente. La + generatrice `θ = 0` è quella del sensore; l'offset circonferenziale delle sorgenti + (`offset_y_percorso_m`, misurato lungo l'arco) è riferito ad essa. +- Lo **spessore non è discretizzato**: la parete è termicamente sottile, quindi la + faccia esterna e quella interna hanno la stessa temperatura. 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). +interna nel punto `x` configurato (default 50 mm, al centro della lunghezza) sulla +generatrice `θ = 0`. ## Modello fisico -Non è una simulazione FEM elettromagnetica + termica completa: è un generatore -pratico di dataset. La catena di approssimazioni è la seguente. +Il campo elettromagnetico non è simulato: la parte termica è però una FEM vera. La +catena di modellazione è la seguente. -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. +1. **Shell termica** — per una parete sottile con temperatura uniforme nello spessore + l'equazione del calore sulla superficie cilindrica è -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. + ```text + ρ cp s ∂T/∂t = k s ( ∂²T/∂x² + (1/R²) ∂²T/∂θ² ) + q'' − q_conv + ``` -3. **Diffusione 2D del calore** — l'equazione del calore è integrata nella sezione - `(x, z)` con volumi finiti ed Eulero implicito (incondizionatamente stabile). + La conduzione **assiale** e quella **circonferenziale** sono entrambe risolte. + Discretizzando con elementi quadrilateri a 4 nodi e una temperatura per nodo si + ottiene `C Ṫ + K T = f(t)`, con -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. + ```text + Cₑ = ∫ ρ cp s Nᵀ N dA (matrice consistente) + Kₑ_cond = ∫ Bᵀ D B s dA D = k·I (materiale isotropo) + ``` + + La conduzione tra elementi adiacenti nasce dall'assemblaggio: due elementi + confinanti condividono i nodi del bordo comune, quindi condividono gli stessi + gradi di libertà termici. Non serve alcuna legge di scambio tra elementi dello + stesso materiale. + +2. **Convezione** — nelle shell la faccia esterna e quella interna sono **entrambe** + superfici fisiche esposte, quindi ogni elemento scambia su tutta la sua area: + + ```text + Kₑ_conv = ∫ (h_est + h_int) Nᵀ N dA fₑ_conv = ∫ (h_est + h_int) T∞ Nᵀ dA + ``` + + I due bordi assiali espongono solo lo spessore della lamiera: + `Kₑ_bordo = ∫ h_b s Nᵀ N dl`. Con `s = 0.18 mm` su `100 mm` di lunghezza il + contributo dei bordi vale ~0.2% dell'area ed è quindi marginale. + +3. **Sorgenti gaussiane in moto** — un gruppo di `numero_sorgenti` sorgenti + equidistanti (`distanza_sorgenti_m`) si muove rigidamente a velocità costante lungo + la generatrice individuata da `offset_y_percorso_m`. Il flusso sulla superficie + esterna è + + ```text + q''(x, θ, t) = q_picco · efficienza · Σᵢ exp( −((x − xᵢ(t))² + Δy²) / (2σ²) ) + Δy = R · wrap(θ − θ_s) + ``` + + La distanza circonferenziale è **riavvolta** su `2πR`: la sorgente scalda anche i + punti raggiunti dall'altro verso della circonferenza. Il carico nodale è + `fₑ_src = ∫ Nᵀ q'' dA`, integrato con la stessa quadratura di Gauss delle matrici: + la posizione circonferenziale dell'induttore è quindi simulata davvero, non + rappresentata da un coefficiente di attenuazione. + +4. **Skin depth** — `δ = √(2ρₑ/(ωμ))`, calcolata dalla resistività elettrica e dalla + permeabilità del materiale, oppure imposta con `skin_depth_fissa_m`. Poiché la + shell non discretizza lo spessore, la skin depth **non entra nella soluzione**: è + registrata in `metadata.csv` come diagnostica. Per la banda stagnata a 20 kHz vale + 0.142 mm, confrontabile con lo spessore di 0.18 mm — cioè la deposizione è già + quasi uniforme attraverso la parete. 5. **Temperatura iniziale** — il campo parte uniformemente alla temperatura ambiente del run (che è randomizzata, quindi varia run per run). -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`). +6. **Sensore realistico** — la temperatura nel punto osservato è interpolata con le + funzioni di forma dell'elemento che lo contiene; la lettura vi aggiunge inerzia del + primo ordine (`costante_tempo_s`), rumore gaussiano (`rumore_std_C`) e + quantizzazione (`quantizzazione_C`). + +### Perché una sola temperatura nello spessore + +La parete è termicamente sottile con ampio margine. Il numero di Fourier che confronta +il tempo di diffusione attraverso lo spessore con il tempo di transito della sorgente, + +```text +Fo = α · (σ/|v|) / s² +``` + +vale **5.4** con i parametri di default (ed è riportato per ogni run in +`metadata.csv`). Il modello 2D `(x, z)` precedente, che risolveva lo spessore con 15 +nodi, dava un salto esterno–interno di **0.002 °C al picco** su una sovratemperatura +di ~200 °C, e 0.53 °C come massimo su tutto il run. Risolvere lo spessore non aggiunge +informazione a questi parametri. + +Il criterio per sapere quando servirebbe una formulazione multistrato è `Fo ≈ 1`, cioè +uno spessore critico `s* = √(α σ/|v|)`: + +| materiale | Fo a 0.18 mm | spessore critico | +|--------------------------|--------------|------------------| +| alluminio | 31.3 | 1.01 mm | +| ferro | 8.4 | 0.52 mm | +| banda stagnata / acciaio dolce | 5.4 | 0.42 mm | +| inox 304 | 1.5 | 0.22 mm | ## 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 strutturata `n_elementi_x × n_elementi_theta` (default 100 × 110, cioè elementi + di 1.0 × 2.0 mm) per **11110 nodi**; la periodicità circonferenziale è strutturale, + non un vincolo imposto: non esiste un nodo separato per `θ = 2π`. +- Quadratura di Gauss 2×2 per tutte le matrici elementari e per il vettore della + sorgente. La mesh è uniforme, quindi tutti gli elementi sono lo stesso rettangolo: + le matrici elementari si calcolano una volta e si replicano. +- Eulero implicito (incondizionatamente stabile) con passo `dt_interno_s` + (default 1 ms), più fine del periodo di campionamento CSV. +- Tutti i termini sono lineari e costanti nel run: `A = C/dt + K_cond + K_conv` viene + **fattorizzata LU una sola volta per run** (`scipy.sparse.linalg.splu`); ogni passo + risolve solo i sistemi triangolari. `A` è simmetrica definita positiva, quindi si usa + l'ordinamento `MMD_AT_PLUS_A`, che dimezza circa il riempimento rispetto al COLAMD di + default. Un run da 30 s simulati richiede ~65 s di calcolo; i run sono generati in + parallelo su più processi. - Se in futuro si introducessero proprietà dipendenti dalla temperatura, la matrice andrebbe ricostruita e rifattorizzata a ogni passo. +### Validazione + +`python valida_shell.py` verifica il solutore su casi con soluzione nota. Va rieseguito +dopo ogni modifica alla fisica o all'assemblaggio. Risultati con la mesh di default: + +| Verifica | Esito | +|-------------------------------------------------|------------------------------------------| +| capacità totale = `ρ cp s A` | esatta a precisione macchina | +| somma righe di `K_cond` = 0 (T costante → flusso nullo) | 4·10⁻¹⁸ | +| simmetria e semidefinita positività di `K_cond` | verificate | +| raffreddamento convettivo uniforme vs `exp(−t/τ)` | errore 0.011 °C, pari a `dt/(2τ)` | +| conservazione dell'energia (adiabatico) | errore relativo 4·10⁻¹⁴ | +| `Σ f_src` vs potenza analitica `2πσ² q_max` | rapporto 1.0000; **0.5000** con sorgente centrata sul bordo | +| modi sinusoidali circonferenziali `n = 1, 2, 4` vs `exp(−α n²t/R²)` | errore 0.04% – 0.17% | +| convergenza al raffinamento (modo `n = 4`) | 2.75% → 1.07% → 0.17% → 0.05% da 30 a 220 elementi | +| simmetria della diffusione attorno a `θ = 0` | 2·10⁻¹⁴ | + +L'ultima riga è la prova che la periodicità funziona davvero: un impulso su `θ = 0` +diffonde in modo identico nei due versi solo se la connettività riavvolge la +circonferenza. + ## 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) +config.py tutti i parametri di simulazione +materials.py proprietà termofisiche ed elettriche dei materiali +mesh.py nodi, connettività Q4, periodicità circonferenziale +elementi_shell.py funzioni di forma, quadratura, matrici del singolo elemento +assemblaggio.py assemblaggio delle matrici globali e del vettore sorgente +sorgente.py moto degli induttori e distribuzione gaussiana del flusso +solutore.py integrazione implicita e fattorizzazione sparsa +sensore.py interpolazione, inerzia, rumore, quantizzazione +simulate.py orchestrazione dei run + generazione dei CSV +valida_shell.py verifiche del solutore su casi con soluzione nota +plot_csv.py grafici rapidi (temperature e flusso) del primo run +plot_animazione.py animazione della superficie sviluppata T(x, θ) +plot_animazione_3d.py animazione 3D isometrica del cilindro +dataset/ output generato da simulate.py (ricreato a ogni esecuzione) ``` ## Installazione @@ -127,19 +213,50 @@ 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 +# animazione della superficie sviluppata durante il passaggio delle sorgenti python plot_animazione.py + +# animazione 3D isometrica del cilindro +python plot_animazione_3d.py ``` Gli script di visualizzazione aprono finestre interattive (backend Qt); se il backend non è interattivo (es. sessione senza display) salvano automaticamente PNG/GIF in `dataset/`. -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`. +Entrambe le animazioni riproducono esattamente la fisica di `run_0001` (stesso seed). +Quella 2D mostra tre pannelli allineati: il profilo di flusso `q(x)` con le sorgenti in +transito, il campo di temperatura sulla superficie sviluppata in piano (asse verticale += arco circonferenziale, con il sensore a quota 0), e la temperatura nel punto +osservato (vera e con inerzia del sensore). Quella 3D disegna il campo calcolato sulla +superficie del cilindro — non è più una ricostruzione, il modello risolve davvero la +distribuzione circonferenziale. Finestra temporale e cadenza dei fotogrammi si regolano +con le costanti in testa a `plot_animazione.py`. + +La mappa sviluppata usa `pcolormesh` con shading Gouraud, che interpola linearmente tra +i valori nodali: è esattamente ciò che fanno le funzioni di forma bilineari del Q4, +quindi si guarda la soluzione dell'elemento finito e non un'interpolazione grafica. + +### Vedere gli elementi shell + +Entrambe le animazioni hanno una costante `MOSTRA_MESH` che disegna i bordi degli +elementi. È `False` di default perché con la mesh di produzione gli elementi sono +1 × 2 mm, più fitti dei pixel disponibili, e il reticolo coprirebbe il campo. Per +ispezionare la mesh servono due modifiche insieme: + +```python +# in plot_animazione.py e/o plot_animazione_3d.py +MOSTRA_MESH = True + +# in config.py, dentro FASCETTA +"n_elementi_x": 25, +"n_elementi_theta": 28, +``` + +Così si vedono i singoli quadrilateri, in piano e avvolti sul cilindro. Attenzione: con +28 elementi sulla circonferenza l'arco vale 7.9 mm contro un `sigma_punto_m` di 12 mm, +cioè meno di due elementi per sigma — la mesh grossolana serve a **guardare la +discretizzazione**, non a produrre risultati (vedi la tabella di convergenza sopra). ## Configurazione @@ -148,8 +265,8 @@ 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 | +| `FASCETTA` | diametro medio, lunghezza, spessore, densità della mesh, materiale | +| `ARIA` | temperatura ambiente e coefficienti di convezione di 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 run | @@ -191,26 +308,30 @@ Punti da conoscere: | `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_max_fascetta_C` | temperatura massima su tutta la superficie in quell'istante | +| `T_media_fascetta_C` | temperatura media su tutta la superficie in quell'istante | | `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. +Contiene tutti i parametri effettivi del run (geometria, dimensioni della mesh, +coefficienti di scambio, parametri delle sorgenti e del sensore, valori randomizzati), +le diagnostiche `skin_depth_m` e `numero_fourier_spessore`, e le temperature di picco: +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). + un'impronta gaussiana moltiplicata per un'efficienza costante. +2. Lo spessore non è discretizzato: il modello non può rappresentare un gradiente tra + faccia esterna e faccia interna. Con i parametri di default il gradiente reale è di + pochi millesimi di grado (vedi `numero_fourier_spessore`), ma con pareti più spesse + di ~0.4 mm, materiali a bassa diffusività o sorgenti molto più veloci servirebbe una + shell multistrato. 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 +4. Irraggiamento non modellato: a ~210 °C le perdite radiative non sono del tutto trascurabili rispetto alla convezione. diff --git a/assemblaggio.py b/assemblaggio.py new file mode 100644 index 0000000..0d74af9 --- /dev/null +++ b/assemblaggio.py @@ -0,0 +1,151 @@ +# Assemblaggio delle matrici globali della shell termica. +# +# Tutti gli elementi della mesh cilindrica sono lo stesso rettangolo, quindi +# ogni matrice elementare si calcola una volta e si replica su tutti gli +# elementi. L'assemblaggio è la somma dei contributi sui gradi di libertà +# condivisi: due elementi confinanti condividono i due nodi del bordo comune, +# ed è questo che genera la conduzione tra elementi adiacenti. Non serve +# nessun accoppiamento aggiuntivo tra elementi dello stesso materiale. + +import numpy as np +import scipy.sparse as sp + +import elementi_shell as el + + +def _assembla_matrice(matrice_elemento: np.ndarray, connettivita: np.ndarray, n_nodi: int): + # Replica una matrice elementare costante su tutti gli elementi e somma i + # contributi sui gradi di libertà condivisi. + n_elementi, n_locali = connettivita.shape + righe = np.repeat(connettivita, n_locali, axis=1).ravel() + colonne = np.tile(connettivita, (1, n_locali)).ravel() + dati = np.tile(matrice_elemento.ravel(), n_elementi) + return sp.coo_matrix( + (dati, (righe, colonne)), shape=(n_nodi, n_nodi) + ).tocsr() + + +def _assembla_vettore(vettore_elemento: np.ndarray, connettivita: np.ndarray, n_nodi: int): + dati = np.tile(vettore_elemento, connettivita.shape[0]) + return np.bincount(connettivita.ravel(), weights=dati, minlength=n_nodi) + + +def assembla_capacita(mesh: dict, materiale: dict, spessore_m: float): + # C = integrale di rho * cp * spessore * N^T N dA. + # + # Matrice consistente: con integrazione temporale implicita conserva + # meglio l'energia rispetto alla versione concentrata sulla diagonale. + capacita_superficiale = ( + materiale["densita_kg_m3"] * materiale["calore_specifico_J_kgK"] * spessore_m + ) + M_e = el.matrice_massa_elemento(mesh["lato_x_m"], mesh["lato_arco_m"]) + return _assembla_matrice( + capacita_superficiale * M_e, mesh["connettivita"], mesh["n_nodi"] + ) + + +def assembla_conduzione(mesh: dict, materiale: dict, spessore_m: float): + # K_cond = integrale di B^T D B * spessore dA, con D isotropa. + K_e = el.matrice_conduzione_elemento( + mesh["lato_x_m"], + mesh["lato_arco_m"], + materiale["conducibilita_termica_W_mK"], + spessore_m, + ) + return _assembla_matrice(K_e, mesh["connettivita"], mesh["n_nodi"]) + + +def assembla_convezione(mesh: dict, aria: dict, spessore_m: float): + # Convezione sulle facce e sui bordi assiali. + # + # Nelle shell la faccia esterna e quella interna sono entrambe superfici + # fisiche esposte all'aria, quindi ogni elemento scambia su tutta la sua + # area con il coefficiente combinato h_esterno + h_interno. I bordi + # assiali (x = 0 e x = lunghezza) espongono invece solo lo spessore della + # lamiera: il loro contributo è proporzionale a h_bordi * spessore ed è + # marginale rispetto a quello delle facce. + T_ambiente = aria["temperatura_ambiente_C"] + h_facce = aria["h_esterno_W_m2K"] + aria["h_interno_W_m2K"] + + M_e = el.matrice_massa_elemento(mesh["lato_x_m"], mesh["lato_arco_m"]) + f_e = el.vettore_carico_uniforme_elemento(mesh["lato_x_m"], mesh["lato_arco_m"]) + + K_conv = _assembla_matrice(h_facce * M_e, mesh["connettivita"], mesh["n_nodi"]) + f_ambiente = _assembla_vettore( + h_facce * T_ambiente * f_e, mesh["connettivita"], mesh["n_nodi"] + ) + + h_bordi_efficace = aria["h_bordi_W_m2K"] * spessore_m + M_bordo = el.matrice_massa_segmento(mesh["lato_arco_m"]) + f_bordo = el.vettore_carico_segmento(mesh["lato_arco_m"]) + + K_conv = K_conv + _assembla_matrice( + h_bordi_efficace * M_bordo, mesh["segmenti_bordo"], mesh["n_nodi"] + ) + f_ambiente = f_ambiente + _assembla_vettore( + h_bordi_efficace * T_ambiente * f_bordo, + mesh["segmenti_bordo"], + mesh["n_nodi"], + ) + + return K_conv, f_ambiente + + +def prepara_assemblatore_sorgente(mesh: dict) -> dict: + # Precalcola quanto serve per assemblare il vettore della sorgente a ogni + # passo temporale: le funzioni di forma nei punti di Gauss pesate per il + # peso di quadratura e lo jacobiano, e le coordinate dei punti di Gauss + # separate per direzione. + det_j = mesh["lato_x_m"] * mesh["lato_arco_m"] / 4.0 + pesi = el.pesi_ai_punti_gauss() + N_gauss = el.matrice_forme_ai_punti_gauss() + offset_x, offset_arco = el.coordinate_locali_punti_gauss( + mesh["lato_x_m"], mesh["lato_arco_m"] + ) + + # Coordinate assolute dei punti di Gauss, per direzione: (n_elementi_x, 2) + # e (n_elementi_theta, 2). + x_gauss = mesh["x_nodi_m"][: mesh["n_elementi_x"], None] + offset_x[None, :] + arco_gauss = mesh["arco_nodi_m"][:, None] + offset_arco[None, :] + + return { + "peso_forme": (pesi * det_j)[:, None] * N_gauss, + "connettivita": mesh["connettivita"], + "n_nodi": mesh["n_nodi"], + "n_elementi_x": mesh["n_elementi_x"], + "n_elementi_theta": mesh["n_elementi_theta"], + "x_gauss_m": x_gauss, + "arco_gauss_m": arco_gauss, + } + + +def assembla_sorgente( + assemblatore: dict, flusso_x_W_m2: np.ndarray, fattore_arco: np.ndarray +) -> np.ndarray: + # f_src = integrale di N^T q'' dA, con q'' fattorizzato come + # q''(x, arco) = flusso_x(x) * fattore_arco(arco). + # + # La fattorizzazione non è un'approssimazione: la gaussiana della + # sorgente è separabile in x e in arco, e tutte le sorgenti del gruppo + # condividono la stessa posizione circonferenziale, quindi il fattore + # circonferenziale è comune e costante nel tempo. + # + # flusso_x_W_m2 ha forma (n_elementi_x, 2) e fattore_arco + # (n_elementi_theta, 2): sono i valori nei punti di Gauss di ciascuna + # direzione. + # + # Per gran parte di un run le sorgenti sono fuori dalla fascetta e il + # fattore assiale è identicamente nullo: in quel caso il vettore è nullo + # e non serve percorrere la mesh. + if not flusso_x_W_m2.any(): + return np.zeros(assemblatore["n_nodi"]) + + q_gauss = ( + flusso_x_W_m2[:, None, :, None] * fattore_arco[None, :, None, :] + ).reshape(-1, 4) + contributi = q_gauss @ assemblatore["peso_forme"] + return np.bincount( + assemblatore["connettivita"].ravel(), + weights=contributi.ravel(), + minlength=assemblatore["n_nodi"], + ) diff --git a/config.py b/config.py index 31c1a0f..b0211d1 100644 --- a/config.py +++ b/config.py @@ -1,26 +1,25 @@ -# Configurazione per il simulatore termico 2D della sezione di una fascetta. +# Configurazione per il simulatore termico a elementi shell di una 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. +# - La fascetta è un anello cilindrico con diametro medio "diametro_medio_mm", +# spessore "spessore_mm" e lunghezza "lunghezza_mm". +# - Il dominio simulato è l'intera superficie cilindrica, discretizzata con +# elementi shell quadrilateri a 4 nodi e una temperatura per nodo. +# - Sistema di coordinate: x = direzione della lunghezza (da 0 a lunghezza), +# theta = coordinata circonferenziale, periodica (theta = 0 coincide con +# theta = 2*pi). Il sensore sta per convenzione su theta = 0. +# - Lo spessore non è discretizzato: la parete è termicamente sottile, quindi +# la faccia esterna e quella interna hanno la stessa temperatura. La +# validità dell'ipotesi è riportata a ogni run come numero di Fourier sullo +# spessore in metadata.csv. +# - Le sorgenti a induzione si muovono in direzione -x sulla superficie +# esterna, lungo la generatrice individuata da "offset_y_percorso_m". +# - La conduzione assiale e quella circonferenziale sono entrambe risolte e +# nascono dall'assemblaggio: elementi adiacenti condividono i nodi del +# bordo comune. +# - Faccia esterna e faccia interna scambiano entrambe per convezione con +# l'aria su tutta l'area; i due bordi assiali scambiano attraverso lo +# spessore della lamiera. # - 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 @@ -46,7 +45,7 @@ SIMULAZIONE = { # Passo di integrazione numerica interna. # Può essere inferiore al periodo di campionamento CSV. - "dt_interno_s": 0.0002, + "dt_interno_s": 0.001, # Seed per la riproducibilità. "seed": 42, @@ -60,8 +59,9 @@ SIMULAZIONE = { } FASCETTA = { - # Diametro della fascetta [mm]. - "diametro_mm": 70.0, + # Diametro medio della fascetta [mm], cioè il diametro della superficie + # su cui giace la mesh shell. + "diametro_medio_mm": 70.0, # Lunghezza della fascetta lungo x [mm]. "lunghezza_mm": 100.0, @@ -69,12 +69,12 @@ FASCETTA = { # Spessore della parete [mm]. "spessore_mm": 0.18, - # Numero di celle del volume finito lungo x (lunghezza). - "n_nodi_x": 100, + # Numero di elementi shell lungo x (lunghezza). + "n_elementi_x": 100, - # Numero di celle del volume finito lungo z (spessore). - # Più nodi = maggiore risoluzione spaziale, simulazione più lenta. - "n_nodi_z": 15, + # Numero di elementi shell lungo la circonferenza. + # Più elementi = maggiore risoluzione spaziale, simulazione più lenta. + "n_elementi_theta": 110, # Deve corrispondere a una chiave in materials.py. "materiale": "banda_stagnata", @@ -84,13 +84,14 @@ 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 assiali (x = 0 e x = lunghezza), + # che espongono all'aria solo lo spessore della lamiera. "h_bordi_W_m2K": 10.0, } @@ -103,9 +104,10 @@ SORGENTE = { # gruppo più lontana dal sensore (quella che lo supera per ultima). "x_fine_m": 5.0, - # Offset circonferenziale (y) tra il percorso delle sorgenti e il punto - # osservato dal sensore. Se diverso da zero, la sorgente passa a lato, - # riducendo il picco. + # Offset circonferenziale, misurato lungo l'arco, tra il percorso delle + # sorgenti e la generatrice osservata dal sensore. Se diverso da zero, la + # sorgente passa a lato e il calore deve diffondere lungo la + # circonferenza per raggiungere il sensore. "offset_y_percorso_m": 0.0, # Velocità delle sorgenti lungo x. Il segno indica il verso di marcia. @@ -128,6 +130,9 @@ SORGENTE = { "efficienza_riscaldamento": 0.35, # Frequenza di induzione usata per stimare la skin depth se skin_depth_fissa_m è None. + # La shell non discretizza lo spessore, quindi la skin depth non entra + # nella soluzione: è calcolata e registrata come diagnostica, per + # verificare che resti confrontabile con lo spessore della lamiera. "frequenza_hz": 20000.0, # Override della skin depth. Usare None per calcolarla dalle proprietà elettriche del materiale. @@ -140,9 +145,11 @@ SORGENTE = { SENSORE = { # Coordinata x del punto della superficie interna osservato dal sensore [mm]. + # La posizione circonferenziale è theta = 0 per convenzione: è l'offset + # delle sorgenti a essere misurato rispetto al sensore. "x_mm": 50.0, - # Distanza del sensore dalla parete interna lungo z [mm]. + # Distanza del sensore dalla parete interna [mm]. # Il sensore è a infrarossi: la distanza non influenza la misura, # è registrata solo come geometria del setup. "distanza_parete_mm": 10.0, diff --git a/elementi_shell.py b/elementi_shell.py new file mode 100644 index 0000000..e0307fb --- /dev/null +++ b/elementi_shell.py @@ -0,0 +1,135 @@ +# Elemento shell termico quadrilatero a 4 nodi (Q4) con una temperatura per +# nodo, cioè temperatura uniforme nello spessore. +# +# L'elemento vive nel piano (x, arco) della superficie sviluppata. Poiché la +# mesh cilindrica è strutturata e uniforme, tutti gli elementi sono lo stesso +# rettangolo lato_x × lato_arco: le matrici elementari si calcolano una volta +# sola e si replicano in fase di assemblaggio. +# +# Le funzioni di forma bilineari sono +# +# N_a(xi, eta) = (1 + xi_a * xi) * (1 + eta_a * eta) / 4 +# +# con (xi_a, eta_a) i vertici dell'elemento di riferimento. La mappa verso +# l'elemento fisico è affine, quindi lo jacobiano è costante: +# det(J) = lato_x * lato_arco / 4. +# +# Integrazione con quadratura di Gauss 2x2, esatta per i polinomi che +# compaiono nelle matrici di massa e conduzione dell'elemento Q4. + +import numpy as np + + +# Vertici dell'elemento di riferimento, nello stesso ordine dei nodi locali +# usati dalla connettività in mesh.py. +NODI_RIFERIMENTO = np.array([(-1.0, -1.0), (1.0, -1.0), (1.0, 1.0), (-1.0, 1.0)]) + +# Quadratura di Gauss a 2 punti per direzione. +PUNTI_GAUSS_1D = np.array([-1.0 / np.sqrt(3.0), 1.0 / np.sqrt(3.0)]) +PESI_GAUSS_1D = np.array([1.0, 1.0]) + + +def funzioni_forma(xi: float, eta: float) -> np.ndarray: + xi_a, eta_a = NODI_RIFERIMENTO[:, 0], NODI_RIFERIMENTO[:, 1] + return 0.25 * (1.0 + xi_a * xi) * (1.0 + eta_a * eta) + + +def derivate_funzioni_forma(xi: float, eta: float) -> tuple[np.ndarray, np.ndarray]: + xi_a, eta_a = NODI_RIFERIMENTO[:, 0], NODI_RIFERIMENTO[:, 1] + dN_dxi = 0.25 * xi_a * (1.0 + eta_a * eta) + dN_deta = 0.25 * eta_a * (1.0 + xi_a * xi) + return dN_dxi, dN_deta + + +def _punti_quadratura_2d(): + # Punti (xi, eta) e pesi della quadratura 2x2, nell'ordine + # g = p * 2 + q con xi = punto p (direzione x) ed eta = punto q + # (direzione circonferenziale). Lo stesso ordine è atteso dai valori del + # flusso passati a vettore_sorgente_elementi(). + for p, xi in enumerate(PUNTI_GAUSS_1D): + for q, eta in enumerate(PUNTI_GAUSS_1D): + yield xi, eta, PESI_GAUSS_1D[p] * PESI_GAUSS_1D[q] + + +def matrice_forme_ai_punti_gauss() -> np.ndarray: + # Matrice (4 punti di Gauss, 4 nodi) con le funzioni di forma valutate + # nei punti di quadratura. + return np.array([funzioni_forma(xi, eta) for xi, eta, _ in _punti_quadratura_2d()]) + + +def pesi_ai_punti_gauss() -> np.ndarray: + # Pesi di quadratura nei 4 punti, nello stesso ordine. + return np.array([peso for _, _, peso in _punti_quadratura_2d()]) + + +def coordinate_locali_punti_gauss(lato_x_m: float, lato_arco_m: float): + # Offset dei punti di Gauss rispetto al vertice di riferimento + # dell'elemento, separati per direzione: 2 valori lungo x e 2 lungo + # l'arco. La separazione è possibile perché l'elemento è un rettangolo. + offset_x = lato_x_m * (1.0 + PUNTI_GAUSS_1D) / 2.0 + offset_arco = lato_arco_m * (1.0 + PUNTI_GAUSS_1D) / 2.0 + return offset_x, offset_arco + + +def matrice_massa_elemento(lato_x_m: float, lato_arco_m: float) -> np.ndarray: + # Integrale di N^T N sull'area dell'elemento. Matrice consistente (non + # concentrata): è la base sia della capacità termica sia della convezione + # sulle facce. + det_j = lato_x_m * lato_arco_m / 4.0 + M = np.zeros((4, 4)) + for xi, eta, peso in _punti_quadratura_2d(): + N = funzioni_forma(xi, eta) + M += peso * det_j * np.outer(N, N) + return M + + +def matrice_conduzione_elemento( + lato_x_m: float, lato_arco_m: float, conducibilita_W_mK: float, spessore_m: float +) -> np.ndarray: + # Integrale di B^T D B * spessore sull'area, con D = k * identità + # (materiale isotropo). Rappresenta la conduzione tangenziale, cioè sia + # quella assiale sia quella circonferenziale: sul cilindro sviluppato le + # due direzioni sono ortogonali e la derivata circonferenziale rispetto + # all'arco equivale a (1/R) d/dtheta. + det_j = lato_x_m * lato_arco_m / 4.0 + K = np.zeros((4, 4)) + for xi, eta, peso in _punti_quadratura_2d(): + dN_dxi, dN_deta = derivate_funzioni_forma(xi, eta) + # Jacobiano diagonale e costante: la derivata fisica è quella + # naturale riscalata dal semilato dell'elemento. + B = np.stack([dN_dxi * (2.0 / lato_x_m), dN_deta * (2.0 / lato_arco_m)]) + K += peso * det_j * conducibilita_W_mK * spessore_m * (B.T @ B) + return K + + +def vettore_carico_uniforme_elemento(lato_x_m: float, lato_arco_m: float) -> np.ndarray: + # Integrale di N^T sull'area: distribuisce ai nodi un carico superficiale + # uniforme unitario. + det_j = lato_x_m * lato_arco_m / 4.0 + f = np.zeros(4) + for xi, eta, peso in _punti_quadratura_2d(): + f += peso * det_j * funzioni_forma(xi, eta) + return f + + +def _funzioni_forma_segmento(xi: float) -> np.ndarray: + return np.array([0.5 * (1.0 - xi), 0.5 * (1.0 + xi)]) + + +def matrice_massa_segmento(lunghezza_m: float) -> np.ndarray: + # Integrale di N^T N su un segmento a 2 nodi, per i bordi assiali. + det_j = lunghezza_m / 2.0 + M = np.zeros((2, 2)) + for xi, peso in zip(PUNTI_GAUSS_1D, PESI_GAUSS_1D): + N = _funzioni_forma_segmento(xi) + M += peso * det_j * np.outer(N, N) + return M + + +def vettore_carico_segmento(lunghezza_m: float) -> np.ndarray: + # Integrale di N^T su un segmento a 2 nodi. + det_j = lunghezza_m / 2.0 + f = np.zeros(2) + for xi, peso in zip(PUNTI_GAUSS_1D, PESI_GAUSS_1D): + f += peso * det_j * _funzioni_forma_segmento(xi) + return f diff --git a/mesh.py b/mesh.py new file mode 100644 index 0000000..0c25fa2 --- /dev/null +++ b/mesh.py @@ -0,0 +1,119 @@ +# Generazione della mesh a elementi shell quadrilateri sulla superficie +# cilindrica della fascetta. +# +# La superficie è parametrizzata da (x, theta): x è la coordinata assiale, +# theta quella circonferenziale. Al posto di theta si usa quasi ovunque +# l'ascissa curvilinea arco = raggio_medio * theta, così la shell diventa un +# rettangolo lunghezza × circonferenza e gli elementi sono tutti identici. +# +# I nodi formano una griglia (n_elementi_x + 1) × n_elementi_theta: in +# direzione circonferenziale non c'è un nodo finale distinto perché theta = 0 +# e theta = 2*pi sono lo stesso nodo. La periodicità è quindi strutturale, +# nasce dalla connettività e non da vincoli imposti a posteriori. +# +# Ordine dei nodi locali di ogni elemento (antiorario nel piano x-arco): +# +# 4 ---- 3 arco +# | | ^ +# | | | +# 1 ---- 2 +---> x + +import numpy as np + + +def costruisci_mesh_cilindrica( + lunghezza_m: float, + raggio_medio_m: float, + n_elementi_x: int, + n_elementi_theta: int, +) -> dict: + # Costruisce nodi e connettività della shell cilindrica. + n_nodi_x = n_elementi_x + 1 + n_nodi = n_nodi_x * n_elementi_theta + n_elementi = n_elementi_x * n_elementi_theta + + lato_x_m = lunghezza_m / n_elementi_x + circonferenza_m = 2.0 * np.pi * raggio_medio_m + lato_arco_m = circonferenza_m / n_elementi_theta + + x_nodi_m = np.arange(n_nodi_x) * lato_x_m + theta_nodi_rad = np.arange(n_elementi_theta) * (2.0 * np.pi / n_elementi_theta) + arco_nodi_m = theta_nodi_rad * raggio_medio_m + + # Connettività: elemento (i, j) collega i nodi (i, j), (i+1, j), + # (i+1, j+1), (i, j+1), con j+1 riavvolto sulla circonferenza. + i_elem = np.repeat(np.arange(n_elementi_x), n_elementi_theta) + j_elem = np.tile(np.arange(n_elementi_theta), n_elementi_x) + j_succ = (j_elem + 1) % n_elementi_theta + + def indice(i, j): + return i * n_elementi_theta + j + + connettivita = np.stack( + [ + indice(i_elem, j_elem), + indice(i_elem + 1, j_elem), + indice(i_elem + 1, j_succ), + indice(i_elem, j_succ), + ], + axis=1, + ) + + # Segmenti dei due bordi assiali (x = 0 e x = lunghezza): anelli chiusi di + # n_elementi_theta segmenti, usati per la convezione sullo spessore. + j = np.arange(n_elementi_theta) + j_dopo = (j + 1) % n_elementi_theta + segmenti_bordo = np.concatenate( + [ + np.stack([indice(0, j), indice(0, j_dopo)], axis=1), + np.stack( + [indice(n_elementi_x, j), indice(n_elementi_x, j_dopo)], axis=1 + ), + ] + ) + + return { + "n_elementi_x": n_elementi_x, + "n_elementi_theta": n_elementi_theta, + "n_elementi": n_elementi, + "n_nodi_x": n_nodi_x, + "n_nodi": n_nodi, + "lunghezza_m": lunghezza_m, + "raggio_medio_m": raggio_medio_m, + "circonferenza_m": circonferenza_m, + "lato_x_m": lato_x_m, + "lato_arco_m": lato_arco_m, + "x_nodi_m": x_nodi_m, + "theta_nodi_rad": theta_nodi_rad, + "arco_nodi_m": arco_nodi_m, + "connettivita": connettivita, + "segmenti_bordo": segmenti_bordo, + "area_totale_m2": lunghezza_m * circonferenza_m, + } + + +def campo_su_griglia(mesh: dict, T: np.ndarray) -> np.ndarray: + # Rimappa il vettore nodale sulla griglia (n_nodi_x, n_elementi_theta), + # comoda per le mappe sviluppate e per le superfici 3D. + return T.reshape(mesh["n_nodi_x"], mesh["n_elementi_theta"]) + + +def coordinate_3d(mesh: dict, chiudi_circonferenza: bool = True): + # Coordinate cartesiane dei nodi per la vista 3D: x lungo l'asse del + # cilindro, y e z sulla sezione circolare. + # + # Con chiudi_circonferenza si ripete la prima colonna in coda, così la + # superficie disegnata non mostra una fessura in theta = 0. È solo una + # necessità di disegno: il nodo ripetuto non è un grado di libertà. + theta = mesh["theta_nodi_rad"] + if chiudi_circonferenza: + theta = np.append(theta, 2.0 * np.pi) + X, Theta = np.meshgrid(mesh["x_nodi_m"], theta, indexing="ij") + R = mesh["raggio_medio_m"] + return X, R * np.sin(Theta), R * np.cos(Theta) + + +def chiudi_campo(campo: np.ndarray) -> np.ndarray: + # Ripete la prima colonna circonferenziale in coda, in accordo con + # coordinate_3d(chiudi_circonferenza=True). + return np.concatenate([campo, campo[:, :1]], axis=1) diff --git a/plot_animazione.py b/plot_animazione.py index 8b63413..a8f685e 100644 --- a/plot_animazione.py +++ b/plot_animazione.py @@ -1,9 +1,14 @@ -# Animazione della sezione della fascetta durante il passaggio delle sorgenti. +# Animazione della superficie 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 profilo di flusso termico q(x) lungo la generatrice delle sorgenti; +# - il campo di temperatura sulla superficie cilindrica sviluppata in piano, +# cioè T(x, arco) con l'arco centrato sulla generatrice del sensore; # - il sensore infrarosso e la temperatura nel punto osservato. +# +# La mappa sviluppata è periodica in verticale: il bordo superiore e quello +# inferiore sono lo stesso punto della fascetta. import random from pathlib import Path @@ -12,13 +17,17 @@ import matplotlib import matplotlib.pyplot as plt import numpy as np from matplotlib.animation import FuncAnimation, PillowWriter +from matplotlib.collections import LineCollection +import sensore as sens +import sorgente as src from config import SIMULAZIONE +from mesh import campo_su_griglia, chiudi_campo from simulate import ( + campo_iniziale, configurazione_randomizzata, - passo_implicito, + passo_termico, prepara_stato_termico, - profilo_flusso_incidente_W_m2, ) # Istante di inizio dei fotogrammi mostrati (la simulazione parte comunque da 0). @@ -33,27 +42,41 @@ DT_FRAME_S = 0.05 # Millisecondi tra i fotogrammi in riproduzione. INTERVALLO_RIPRODUZIONE_MS = 30 +# Disegna i bordi degli elementi shell sopra il campo di temperatura. +# Con la mesh di default (100 x 110, elementi di 1 x 2 mm) il reticolo è più +# fitto dei pixel disponibili e copre il campo: per ispezionare la mesh +# conviene attivarlo e insieme abbassare n_elementi_x e n_elementi_theta in +# config.py, per esempio a 25 x 28. +MOSTRA_MESH = False + 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. + # + # I campi sono restituiti come griglie (n_nodi_x, n_elementi_theta) con la + # colonna dell'arco ruotata in modo da portare la generatrice del sensore + # (theta = 0) al centro: è solo una scelta di visualizzazione. 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"] + stato = prepara_stato_termico(fascetta, aria, sorgente, sensore) + mesh = stato["mesh"] dt = stato["dt_s"] + interpolatore = stato["interpolatore_sensore"] - x_sensore = sensore["x_mm"] / 1000.0 - i_sensore = min(n_x - 1, max(0, int(x_sensore / stato["dx_m"]))) + # Dopo np.roll(..., rotazione) la colonna theta = 0 finisce all'indice + # "rotazione", quindi l'arco con segno della colonna m vale + # (m - rotazione) * lato_arco. + rotazione = mesh["n_elementi_theta"] // 2 + arco_centrato_m = ( + np.arange(mesh["n_elementi_theta"]) - rotazione + ) * mesh["lato_arco_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) + T = campo_iniziale(stato) + T_letta = sens.temperatura_osservata_C(interpolatore, T) tempi, campi, flussi, x_riferimenti = [], [], [], [] T_vere, T_lette = [], [] @@ -61,18 +84,22 @@ def simula_campi(cfg_run: dict) -> dict: 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 + T, info = passo_termico(stato, T, t) + T_vera = sens.temperatura_osservata_C(interpolatore, T) + T_letta = sens.aggiorna_inerzia_C( + T_letta, T_vera, dt, sensore["costante_tempo_s"] + ) if t + 1e-12 >= prossimo_frame_t: + _, q_nodi = src.profilo_flusso_x_W_m2( + sorgente, stato["x_sensore_m"], t, mesh["x_nodi_m"] + ) 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) + campi.append(np.roll(campo_su_griglia(mesh, T), rotazione, axis=1)) + flussi.append(q_nodi) + x_riferimenti.append(info["x_riferimento_m"]) + T_vere.append(T_vera) + T_lette.append(T_letta) prossimo_frame_t += DT_FRAME_S t += dt @@ -84,11 +111,18 @@ def simula_campi(cfg_run: dict) -> dict: "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"], + "mesh": mesh, + "rotazione": rotazione, + "x_nodi_mm": mesh["x_nodi_m"] * 1000.0, + "arco_centrato_mm": arco_centrato_m * 1000.0, + "arco_chiuso_mm": np.append( + arco_centrato_m, arco_centrato_m[-1] + mesh["lato_arco_m"] + ) * 1000.0, "lunghezza_mm": fascetta["lunghezza_mm"], "x_sensore_mm": sensore["x_mm"], - "sorgente": cfg_run["sorgente"], + "offset_sorgente_mm": sorgente["offset_y_percorso_m"] * 1000.0, + "sorgente": sorgente, + "T_ambiente_C": aria["temperatura_ambiente_C"], } @@ -102,23 +136,27 @@ def main() -> None: 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) + arco_mm = dati["arco_chiuso_mm"] + arco_min_mm, arco_max_mm = arco_mm[0], arco_mm[-1] - q_max_MW = max(q.max() for q in dati["flussi"]) / 1e6 + # Il minimo evita un asse degenere se la finestra dell'animazione si + # chiude prima che le sorgenti entrino nella fascetta. + q_max_MW = max(max(q.max() for q in dati["flussi"]) / 1e6, 1e-3) 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 + offset_mm = dati["offset_sorgente_mm"] - fig, (ax_flusso, ax_sezione, ax_storia) = plt.subplots( + fig, (ax_flusso, ax_superficie, 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") + fig.suptitle("Superficie sviluppata della fascetta: sorgenti in transito e sensore") - # Pannello 1: profilo di flusso sul lato esterno e posizioni delle sorgenti. + # Pannello 1: profilo di flusso lungo x 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 @@ -135,37 +173,50 @@ def main() -> None: 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", + # Pannello 2: temperatura sulla superficie sviluppata. L'asse verticale è + # l'arco circonferenziale, con 0 sulla generatrice del sensore. + # + # pcolormesh con shading "gouraud" interpola linearmente tra i valori + # nodali: è esattamente ciò che fanno le funzioni di forma bilineari + # dell'elemento Q4, quindi il disegno è la soluzione dell'elemento finito + # e non un'interpolazione grafica arbitraria. + Xg, Yg = np.meshgrid(dati["x_nodi_mm"], dati["arco_chiuso_mm"], indexing="ij") + mappa = ax_superficie.pcolormesh( + Xg, Yg, + chiudi_campo(dati["campi"][indice_inizio]), + shading="gouraud", cmap="inferno", - vmin=cfg_run["aria"]["temperatura_ambiente_C"], + vmin=dati["T_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, + if MOSTRA_MESH: + # Bordi degli elementi shell: linee lungo x a ogni nodo + # circonferenziale e viceversa. + segmenti = [ + [(x, arco_min_mm), (x, arco_max_mm)] for x in dati["x_nodi_mm"] + ] + [ + [(0.0, a), (lunghezza_mm, a)] for a in dati["arco_chiuso_mm"] + ] + ax_superficie.add_collection( + LineCollection(segmenti, colors="white", linewidths=0.5, alpha=0.55) + ) + ax_superficie.set_xlim(*x_vista_mm) + ax_superficie.set_ylim(arco_min_mm, arco_max_mm) + ax_superficie.set_ylabel("arco circonferenziale [mm]") + ax_superficie.set_xlabel("x [mm]") + ax_superficie.plot( + [dati["x_sensore_mm"]], [0.0], "^", color="tab:blue", markersize=11 ) - ax_sezione.text( - x_sens + 3, 2.4 * spessore_mm, "sensore IR", - color="tab:blue", va="center", + ax_superficie.text( + dati["x_sensore_mm"] + 3, 0.0, "sensore IR", color="tab:blue", va="center" + ) + ax_superficie.axhline( + offset_mm, color="tab:cyan", linewidth=0.8, linestyle="--", alpha=0.7 ) # La colorbar è agganciata a tutti i pannelli per non restringere solo - # quello della sezione, mantenendo allineati gli assi x. + # quello della superficie, mantenendo allineati gli assi x. barra = fig.colorbar( - immagine, ax=(ax_flusso, ax_sezione, ax_storia), pad=0.02, aspect=35 + mappa, ax=(ax_flusso, ax_superficie, ax_storia), pad=0.02, aspect=35 ) barra.set_label("T [°C]") @@ -186,7 +237,7 @@ def main() -> None: k = indice_inizio + frame t = tempi[k] - linea_flusso.set_data(dati["x_centri_mm"], dati["flussi"][k] / 1e6) + linea_flusso.set_data(dati["x_nodi_mm"], dati["flussi"][k] / 1e6) x_sorgenti_mm = ( dati["x_riferimenti"][k] * 1000.0 @@ -198,7 +249,7 @@ def main() -> None: np.full(int(visibili.sum()), q_max_MW * 1.12), ) - immagine.set_data(dati["campi"][k].T) + mappa.set_array(chiudi_campo(dati["campi"][k])) linea_vera.set_data(tempi[: k + 1], dati["T_vere"][: k + 1]) linea_letta.set_data(tempi[: k + 1], dati["T_lette"][: k + 1]) @@ -206,7 +257,7 @@ def main() -> None: testo_tempo.set_text(f"t = {t:.3f} s") return ( - linea_flusso, marker_sorgenti, immagine, + linea_flusso, marker_sorgenti, mappa, linea_vera, linea_letta, cursore, testo_tempo, ) @@ -216,7 +267,7 @@ def main() -> None: # 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" + percorso = Path("dataset") / "animazione_superficie.gif" animazione.save( percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS) ) diff --git a/plot_animazione_3d.py b/plot_animazione_3d.py index 9207912..082fc69 100644 --- a/plot_animazione_3d.py +++ b/plot_animazione_3d.py @@ -1,26 +1,26 @@ # Animazione 3D del barattolo in vista isometrica: colormap della temperatura -# sulla superficie esterna. +# sulla superficie cilindrica. # -# Il modello risolve solo la sezione x-z (vedi CLAUDE.md): la coordinata -# circonferenziale y non è simulata, è collassata in un'attenuazione -# gaussiana del flusso. Per estrudere il campo attorno alla circonferenza si -# riusa la stessa gaussiana: la sovratemperatura rispetto al piano y=0 (dove -# si trova il sensore) viene scalata in funzione della distanza -# circonferenziale dal percorso delle sorgenti. È quindi una ricostruzione -# visiva, non un calcolo di diffusione in y. +# Il modello shell risolve la temperatura sull'intera superficie, quindi il +# campo disegnato qui è quello calcolato: non c'è più nessuna ricostruzione +# della distribuzione circonferenziale. La faccia esterna e quella interna +# hanno la stessa temperatura (una temperatura per nodo), quindi la superficie +# mostrata vale per entrambe. import random +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 pathlib import Path -from config import FASCETTA, SIMULAZIONE +import sorgente as src +from config import SIMULAZIONE +from mesh import chiudi_campo, coordinate_3d from plot_animazione import DT_FRAME_S as DT_FRAME_S_BASE -from plot_animazione import T_FINE_ANIMAZIONE_S, T_INIZIO_ANIMAZIONE_S, simula_campi +from plot_animazione import T_INIZIO_ANIMAZIONE_S, simula_campi from simulate import configurazione_randomizzata # Tempo simulato tra un fotogramma e il successivo dell'animazione 3D. Più @@ -31,19 +31,12 @@ DT_FRAME_S = 0.15 # Millisecondi tra i fotogrammi in riproduzione. INTERVALLO_RIPRODUZIONE_MS = 60 -# Numero di punti lungo la circonferenza per l'estrusione della superficie. -N_THETA = 72 - - -def attenuazione_circonferenziale( - y_m: np.ndarray, offset_y_m: float, sigma_m: float -) -> np.ndarray: - # Fattore che scala la sovratemperatura (T - T_ambiente) del piano y=0 - # in funzione della distanza circonferenziale y dal piano stesso, - # normalizzato in modo da valere 1 in y=0. - sigma = max(sigma_m, 1e-9) - esponente = -0.5 * ((y_m - offset_y_m) ** 2 - offset_y_m ** 2) / (sigma * sigma) - return np.exp(esponente) +# Disegna i bordi degli elementi shell sulla superficie del cilindro. +# Con la mesh di default (100 x 110, elementi di 1 x 2 mm) il reticolo satura +# il disegno e sbiadisce i colori: per ispezionare la mesh conviene attivarlo +# e insieme abbassare n_elementi_x e n_elementi_theta in config.py, per +# esempio a 25 x 28. +MOSTRA_MESH = False def main() -> None: @@ -56,26 +49,19 @@ def main() -> None: passo = max(1, round(DT_FRAME_S / DT_FRAME_S_BASE)) indici_frame = list(range(indice_inizio, len(tempi), passo)) - T_ambiente = cfg_run["aria"]["temperatura_ambiente_C"] - sorgente = cfg_run["sorgente"] - sigma_m = sorgente["sigma_punto_m"] - offset_y_m = sorgente["offset_y_percorso_m"] - numero_sorgenti = sorgente.get("numero_sorgenti", 1) - distanza_m = sorgente.get("distanza_sorgenti_m", 0.0) + mesh = dati["mesh"] + # simula_campi ruota le colonne per centrare il sensore nella mappa + # sviluppata: qui serve l'ordine originale della mesh. + rotazione = dati["rotazione"] - raggio_m = (FASCETTA["diametro_mm"] / 1000.0) / 2.0 - x_centri_m = dati["x_centri_mm"] / 1000.0 - lunghezza_m = dati["lunghezza_mm"] / 1000.0 + T_ambiente = dati["T_ambiente_C"] + sorgente = dati["sorgente"] + raggio_m = mesh["raggio_medio_m"] + lunghezza_m = mesh["lunghezza_m"] - theta = np.linspace(-np.pi, np.pi, N_THETA) - y_circ_m = theta * raggio_m - attenuazione = attenuazione_circonferenziale(y_circ_m, offset_y_m, sigma_m) + Xm, Ym, Zm = coordinate_3d(mesh) - Xm, Thetam = np.meshgrid(x_centri_m, theta) - Ym = raggio_m * np.sin(Thetam) - Zm = raggio_m * np.cos(Thetam) - - T_max = max(c[:, 0].max() for c in dati["campi"]) + T_max = max(c.max() for c in dati["campi"]) # vmin più basso della temperatura ambiente reale: altrimenti le zone # fredde cadrebbero sul nero puro di "inferno" e, essendo lo shading # moltiplicativo, nessuna illuminazione basterebbe a renderle visibili. @@ -101,25 +87,29 @@ def main() -> None: barra = fig.colorbar(mappabile, ax=ax, shrink=0.6, pad=0.05) barra.set_label("T [°C]") + theta_sorgente = src.theta_percorso_rad(sorgente, raggio_m) + def disegna_frame(k: int): ax.cla() ax.view_init(elev=35.264, azim=45) ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m)) ax.set_axis_off() - T_lato_esterno = dati["campi"][k][:, 0] - T_superficie = T_ambiente + (T_lato_esterno[None, :] - T_ambiente) * attenuazione[:, None] - colori = cmap(norm(T_superficie)) + campo = np.roll(dati["campi"][k], -rotazione, axis=1) + colori = cmap(norm(chiudi_campo(campo))) + # Ogni faccia disegnata è un elemento shell della mesh: con + # MOSTRA_MESH se ne vedono i bordi. + bordi = {"edgecolor": (1.0, 1.0, 1.0, 0.25), "linewidth": 0.15} ax.plot_surface( Xm, Ym, Zm, facecolors=colori, rstride=1, cstride=1, - antialiased=False, shade=True, lightsource=lightsource, linewidth=0, + antialiased=False, shade=True, lightsource=lightsource, + **(bordi if MOSTRA_MESH else {"linewidth": 0}), ) - x_sorgenti_m = dati["x_riferimenti"][k] + np.arange(numero_sorgenti) * distanza_m + x_sorgenti_m = src.posizioni_sorgenti_m(sorgente, dati["x_riferimenti"][k]) visibili = (x_sorgenti_m >= 0.0) & (x_sorgenti_m <= lunghezza_m) if visibili.any(): - theta_sorgente = offset_y_m / raggio_m ax.scatter( x_sorgenti_m[visibili], np.full(int(visibili.sum()), raggio_m * np.sin(theta_sorgente) * 1.05), diff --git a/sensore.py b/sensore.py new file mode 100644 index 0000000..5861355 --- /dev/null +++ b/sensore.py @@ -0,0 +1,60 @@ +# Pirometro a infrarossi montato all'interno della fascetta. +# +# Il sensore osserva senza contatto un punto fisso della superficie interna, +# individuato dalla coordinata assiale x e dalla posizione circonferenziale +# theta = 0, che per convenzione è la generatrice di riferimento del modello. +# +# Con una shell a una temperatura per nodo la faccia interna e quella esterna +# hanno la stessa temperatura: il valore osservato è la temperatura della +# shell nel punto, interpolata con le funzioni di forma dell'elemento che lo +# contiene. L'approssimazione regge finché la parete è termicamente sottile +# rispetto al transito della sorgente (vedi numero_fourier_spessore in +# simulate.py). +# +# Alla temperatura vera si sovrappone la catena di misura: inerzia del primo +# ordine, rumore gaussiano e quantizzazione della lettura. + +import numpy as np + +import elementi_shell as el + + +def prepara_interpolatore(mesh: dict, x_m: float, theta_rad: float) -> dict: + # Individua l'elemento che contiene il punto osservato e i pesi nodali + # con cui ricostruirvi la temperatura. + n_x = mesh["n_elementi_x"] + n_theta = mesh["n_elementi_theta"] + + x_limitato = min(max(x_m, 0.0), mesh["lunghezza_m"]) + i = min(n_x - 1, max(0, int(x_limitato / mesh["lato_x_m"]))) + xi = 2.0 * (x_limitato - mesh["x_nodi_m"][i]) / mesh["lato_x_m"] - 1.0 + + arco = (theta_rad * mesh["raggio_medio_m"]) % mesh["circonferenza_m"] + j = int(arco / mesh["lato_arco_m"]) % n_theta + eta = 2.0 * (arco - mesh["arco_nodi_m"][j]) / mesh["lato_arco_m"] - 1.0 + + elemento = i * n_theta + j + return { + "nodi": mesh["connettivita"][elemento], + "pesi": el.funzioni_forma(xi, eta), + "x_m": x_limitato, + "theta_rad": theta_rad, + } + + +def temperatura_osservata_C(interpolatore: dict, T: np.ndarray) -> float: + return float(T[interpolatore["nodi"]] @ interpolatore["pesi"]) + + +def aggiorna_inerzia_C( + T_letta_C: float, T_vera_C: float, dt_s: float, costante_tempo_s: float +) -> float: + # Risposta del primo ordine del sensore. + tau = max(costante_tempo_s, 1e-9) + return T_letta_C + (T_vera_C - T_letta_C) * dt_s / tau + + +def quantizza(valore: float, passo: float) -> float: + if passo <= 0.0: + return valore + return round(valore / passo) * passo diff --git a/simulate.py b/simulate.py index 7b482b0..de53713 100644 --- a/simulate.py +++ b/simulate.py @@ -1,3 +1,11 @@ +# Motore principale: genera N run randomizzati della fascetta riscaldata da +# sorgenti a induzione in movimento e scrive i CSV del dataset. +# +# Il modello termico è una shell a elementi finiti quadrilateri sull'intera +# superficie cilindrica, con una temperatura per nodo. La pipeline di un run è +# in prepara_stato_termico() e passo_termico(), condivise con i moduli di +# animazione: ogni modifica alla fisica va fatta lì, non duplicata. + import csv import math import os @@ -8,11 +16,14 @@ from copy import deepcopy from pathlib import Path import numpy as np -import scipy.sparse as sp -from scipy.sparse.linalg import splu +import assemblaggio +import sensore as sens +import solutore as sol +import sorgente as src from config import ARIA, FASCETTA, RANDOMIZZAZIONE, SENSORE, SIMULAZIONE, SORGENTE from materials import MATERIALI +from mesh import costruisci_mesh_cilindrica MU0 = 4.0 * math.pi * 1e-7 @@ -31,241 +42,120 @@ def calcola_skin_depth_m(materiale: dict, frequenza_hz: float) -> float: 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. +def numero_fourier_spessore( + materiale: dict, spessore_m: float, sorgente: dict +) -> float: + # Verifica quantitativa dell'ipotesi di parete termicamente sottile: # - # 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]: + # Fo = alpha * (sigma / |v|) / spessore^2 # - # 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) + # confronta il tempo di diffusione attraverso lo spessore con il tempo di + # transito della sorgente. Per Fo >> 1 lo spessore si equilibra molto + # prima che la sorgente sia passata e una shell a una temperatura per + # nodo è adeguata; scendendo verso Fo ~ 1 servirebbe una formulazione + # multistrato. + alpha = materiale["conducibilita_termica_W_mK"] / ( + materiale["densita_kg_m3"] * materiale["calore_specifico_J_kgK"] ) - return splu(sp.csc_matrix(A)) + velocita = abs(sorgente["velocita_m_s"]) + if velocita <= 0.0: + return float("inf") + tempo_transito = sorgente["sigma_punto_m"] / velocita + return alpha * tempo_transito / (spessore_m * spessore_m) -def prepara_stato_termico(fascetta: dict, aria: dict, sorgente: dict) -> dict: - # Prepara griglia, coefficienti e solutore fattorizzato per un run: +def prepara_stato_termico( + fascetta: dict, aria: dict, sorgente: dict, sensore: dict +) -> dict: + # Costruisce mesh, matrici globali e solutore fattorizzato per un run: # tutto ciò che resta costante durante l'integrazione temporale. materiale = MATERIALI[fascetta["materiale"]] - lunghezza = fascetta["lunghezza_mm"] / 1000.0 - spessore = fascetta["spessore_mm"] / 1000.0 - n_x = fascetta["n_nodi_x"] - n_z = fascetta["n_nodi_z"] - dx = lunghezza / n_x - dz = spessore / n_z + lunghezza_m = fascetta["lunghezza_mm"] / 1000.0 + spessore_m = fascetta["spessore_mm"] / 1000.0 + raggio_medio_m = fascetta["diametro_medio_mm"] / 2000.0 - dt = SIMULAZIONE["dt_interno_s"] + mesh = costruisci_mesh_cilindrica( + lunghezza_m=lunghezza_m, + raggio_medio_m=raggio_medio_m, + n_elementi_x=fascetta["n_elementi_x"], + n_elementi_theta=fascetta["n_elementi_theta"], + ) + + C = assemblaggio.assembla_capacita(mesh, materiale, spessore_m) + K_cond = assemblaggio.assembla_conduzione(mesh, materiale, spessore_m) + K_conv, f_ambiente = assemblaggio.assembla_convezione(mesh, aria, spessore_m) + + dt_s = SIMULAZIONE["dt_interno_s"] + solutore = sol.costruisci_solutore(C, K_cond + K_conv, dt_s) + + assemblatore = assemblaggio.prepara_assemblatore_sorgente(mesh) + + # Il fattore circonferenziale del flusso non dipende dal tempo: le + # sorgenti traslano solo lungo x. + fattore_arco = src.fattore_circonferenziale( + sorgente, raggio_medio_m, assemblatore["arco_gauss_m"] + ) + attenuazione_al_sensore = float( + src.fattore_circonferenziale(sorgente, raggio_medio_m, np.zeros(1))[0] + ) if sorgente["skin_depth_fissa_m"] is None: - skin_depth = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"]) + skin_depth_m = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"]) else: - skin_depth = float(sorgente["skin_depth_fissa_m"]) + skin_depth_m = 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 + x_sensore_m = sensore["x_mm"] / 1000.0 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), + "mesh": mesh, + "materiale": materiale, + "spessore_m": spessore_m, + "raggio_medio_m": raggio_medio_m, + "dt_s": dt_s, + "solutore": solutore, + "assemblatore": assemblatore, + "f_ambiente": f_ambiente, + "fattore_arco": fattore_arco, + "attenuazione_al_sensore": attenuazione_al_sensore, + "x_sensore_m": x_sensore_m, + "interpolatore_sensore": sens.prepara_interpolatore(mesh, x_sensore_m, 0.0), + "sorgente": sorgente, "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"], + "skin_depth_m": skin_depth_m, + "numero_fourier_spessore": numero_fourier_spessore( + materiale, spessore_m, sorgente ), } -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, :] +def campo_iniziale(stato: dict) -> np.ndarray: + return np.full(stato["mesh"]["n_nodi"], stato["T_ambiente_C"], dtype=float) + + +def passo_termico(stato: dict, T: np.ndarray, t_s: float) -> tuple[np.ndarray, dict]: + # Avanza il campo di temperatura di un passo dt e restituisce anche lo + # stato della sorgente in quell'istante. + sorgente = stato["sorgente"] + x_riferimento_m, flusso_x = src.profilo_flusso_x_W_m2( + sorgente, stato["x_sensore_m"], t_s, stato["assemblatore"]["x_gauss_m"] + ) + f_sorgente = assemblaggio.assembla_sorgente( + stato["assemblatore"], flusso_x, stato["fattore_arco"] + ) + T_next = sol.passo_implicito( + stato["solutore"], T, stato["f_ambiente"], f_sorgente ) - 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 + _, flusso_al_sensore = src.profilo_flusso_x_W_m2( + sorgente, stato["x_sensore_m"], t_s, np.array([stato["x_sensore_m"]]) + ) + return T_next, { + "x_riferimento_m": x_riferimento_m, + "flusso_al_sensore_W_m2": float( + flusso_al_sensore[0] * stato["attenuazione_al_sensore"] + ), + } def configurazione_randomizzata(indice_run: int, rng: random.Random) -> dict: @@ -323,25 +213,22 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: 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"] + stato = prepara_stato_termico(fascetta, aria, sorgente, sensore) + mesh = stato["mesh"] + interpolatore = stato["interpolatore_sensore"] 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] + T = campo_iniziale(stato) + T_letta = sens.temperatura_osservata_C(interpolatore, T) prossimo_campione_t = 0.0 - T_vera_max = T[i_sensore, -1] - T_misurata_max = T_sensore + T_vera_max = T_letta + T_misurata_max = T_letta + T_fascetta_max = float(T.max()) output_csv.parent.mkdir(parents=True, exist_ok=True) @@ -356,7 +243,8 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: "skin_depth_m", "T_vera_lato_sensore_C", "T_misurata_sensore_C", - "T_lato_caldo_C", + "T_max_fascetta_C", + "T_media_fascetta_C", "T_ambiente_C", "velocita_m_s", "sigma_punto_m", @@ -366,37 +254,41 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: 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) + T, info = passo_termico(stato, T, t) # Temperatura vera della superficie interna nel punto osservato # dal sensore infrarosso. - T_vera_lato_sensore = T[i_sensore, -1] + T_vera = sens.temperatura_osservata_C(interpolatore, T) + T_letta = sens.aggiorna_inerzia_C( + T_letta, T_vera, dt, sensore["costante_tempo_s"] + ) - # 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 + # I massimi delle temperature vere sono aggiornati a ogni passo + # interno, non solo agli istanti di campionamento: servono come + # ground truth e il picco è stretto rispetto al periodo di + # campionamento CSV, che da solo lo taglierebbe di oltre un grado. + T_max = float(T.max()) + T_vera_max = max(T_vera_max, T_vera) + T_fascetta_max = max(T_fascetta_max, T_max) # 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"]) + misurata = T_letta + rng.gauss(0.0, sensore["rumore_std_C"]) + misurata = sens.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"{info['x_riferimento_m']:.9f}", f"{sorgente['offset_y_percorso_m']:.9f}", - f"{q_x[i_sensore]:.6f}", + f"{info['flusso_al_sensore_W_m2']:.6f}", f"{skin_depth:.9e}", - f"{T_vera_lato_sensore:.6f}", + f"{T_vera:.6f}", f"{misurata:.6f}", - f"{T[i_sensore, 0]:.6f}", + f"{T_max:.6f}", + f"{float(T.mean()):.6f}", f"{aria['temperatura_ambiente_C']:.6f}", f"{sorgente['velocita_m_s']:.9f}", f"{sorgente['sigma_punto_m']:.9f}", @@ -411,11 +303,14 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: "id_run": cfg_run["id_run"], "file_csv": str(output_csv.name), "materiale": nome_materiale, - "diametro_m": fascetta["diametro_mm"] / 1000.0, + "diametro_medio_m": fascetta["diametro_medio_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, + "n_elementi_x": mesh["n_elementi_x"], + "n_elementi_theta": mesh["n_elementi_theta"], + "n_nodi": mesh["n_nodi"], + "lato_elemento_x_m": mesh["lato_x_m"], + "lato_elemento_arco_m": mesh["lato_arco_m"], "durata_s": durata, "frequenza_campionamento_hz": SIMULAZIONE["frequenza_campionamento_hz"], "dt_interno_s": dt, @@ -425,9 +320,12 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: "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, + "x_sensore_m": stato["x_sensore_m"], "distanza_sensore_parete_m": sensore["distanza_parete_mm"] / 1000.0, "offset_y_percorso_m": sorgente["offset_y_percorso_m"], + "theta_percorso_rad": src.theta_percorso_rad( + sorgente, stato["raggio_medio_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), @@ -436,11 +334,13 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict: "efficienza_riscaldamento": sorgente["efficienza_riscaldamento"], "frequenza_hz": sorgente["frequenza_hz"], "skin_depth_m": skin_depth, + "numero_fourier_spessore": stato["numero_fourier_spessore"], "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, + "T_max_fascetta_C": T_fascetta_max, } diff --git a/solutore.py b/solutore.py new file mode 100644 index 0000000..a379b26 --- /dev/null +++ b/solutore.py @@ -0,0 +1,39 @@ +# Integrazione temporale implicita del sistema semidiscreto della shell. +# +# Dopo l'assemblaggio il problema è +# +# C * dT/dt + K * T = f(t) con K = K_cond + K_conv +# +# Con Eulero implicito: +# +# (C/dt + K) * T^(n+1) = (C/dt) * T^n + f^(n+1) +# +# La matrice A = C/dt + K dipende solo da mesh, materiale, passo temporale e +# condizioni convettive: nessuno di questi cambia durante un run, quindi si +# fattorizza LU una volta sola e a ogni passo si risolvono solo i sistemi +# triangolari. Se in futuro si introducessero proprietà dipendenti dalla +# temperatura, A andrebbe ricostruita e rifattorizzata a ogni passo. + +import scipy.sparse as sp +from scipy.sparse.linalg import splu + + +def costruisci_solutore(C, K, dt_s: float) -> dict: + # A è simmetrica e definita positiva: capacità, conduzione e convezione + # lo sono tutte. L'ordinamento MMD sul grafo di A + A^T sfrutta la + # simmetria e riduce il riempimento della fattorizzazione rispetto al + # COLAMD di default, pensato per matrici non simmetriche. + A = (C / dt_s + K).tocsc() + return { + "lu": splu(A, permc_spec="MMD_AT_PLUS_A"), + "C_su_dt": (C / dt_s).tocsr(), + "dt_s": dt_s, + } + + +def passo_implicito(solutore: dict, T, f_costante, f_sorgente): + # Avanza il campo di un passo dt. f_costante raccoglie i termini noti + # indipendenti dal tempo (convezione verso l'ambiente su facce e bordi), + # f_sorgente il carico della sorgente all'istante di arrivo. + rhs = solutore["C_su_dt"] @ T + f_costante + f_sorgente + return solutore["lu"].solve(rhs) diff --git a/sorgente.py b/sorgente.py new file mode 100644 index 0000000..20276c9 --- /dev/null +++ b/sorgente.py @@ -0,0 +1,125 @@ +# Sorgenti a induzione in movimento sulla superficie esterna della fascetta. +# +# Il gruppo è formato da più sorgenti equidistanti che traslano insieme lungo +# x con la stessa velocità. Tutte percorrono la stessa generatrice, cioè la +# stessa posizione circonferenziale, definita dall'offset rispetto al punto +# osservato dal sensore. +# +# Il flusso termico incidente su un punto della superficie è +# +# q''(x, arco, t) = q_picco * somma_i exp(-((x - x_i(t))^2) / (2 sigma^2)) +# * exp(-(delta_arco^2) / (2 sigma^2)) +# +# dove delta_arco è la distanza circonferenziale minima tra il punto e il +# percorso delle sorgenti, misurata lungo la circonferenza e riavvolta su +# 2*pi*R: passare da un lato o dall'altro della fascetta è equivalente. +# +# Il flusso è quindi il prodotto di un fattore che dipende solo da x e varia +# nel tempo e di un fattore che dipende solo dall'arco ed è costante: questa +# separazione è sfruttata dall'assemblaggio del vettore sorgente. + +import math + +import numpy as np + + +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 theta_percorso_rad(sorgente: dict, raggio_medio_m: float) -> float: + # Posizione circonferenziale del percorso delle sorgenti, ricavata + # dall'offset in arco rispetto al sensore, che sta per convenzione a + # theta = 0. + return sorgente["offset_y_percorso_m"] / raggio_medio_m + + +def fattore_circonferenziale( + sorgente: dict, raggio_medio_m: float, arco_m: np.ndarray +) -> np.ndarray: + # Attenuazione gaussiana del flusso in funzione della distanza + # circonferenziale dal percorso delle sorgenti. La distanza è riavvolta + # sulla circonferenza, così la sorgente scalda anche i punti raggiunti + # "dall'altra parte" della fascetta. + circonferenza = 2.0 * math.pi * raggio_medio_m + delta = arco_m - sorgente["offset_y_percorso_m"] + delta = (delta + circonferenza / 2.0) % circonferenza - circonferenza / 2.0 + sigma = sorgente["sigma_punto_m"] + return np.exp(-0.5 * (delta * delta) / (sigma * sigma)) + + +def profilo_flusso_x_W_m2( + sorgente: dict, x_sensore_m: float, t_s: float, x_m: np.ndarray +) -> tuple[float, np.ndarray]: + # Restituisce la posizione della sorgente di riferimento e il fattore + # assiale del flusso efficace [W/m²], somma dei contributi di tutte le + # sorgenti attive valutata nei punti x richiesti. + 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"] + ) + + q_x = np.zeros_like(x_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_m - x_i + q_x += q_picco * np.exp(-0.5 * (dx * dx) / (sigma * sigma)) + + return x_riferimento, q_x + + +def posizioni_sorgenti_m(sorgente: dict, x_riferimento_m: float) -> np.ndarray: + # Posizioni lungo x di tutte le sorgenti del gruppo a partire da quella + # di riferimento. + numero_sorgenti = sorgente.get("numero_sorgenti", 1) + distanza = sorgente.get("distanza_sorgenti_m", 0.0) + return x_riferimento_m + np.arange(numero_sorgenti) * distanza diff --git a/valida_shell.py b/valida_shell.py new file mode 100644 index 0000000..87a122a --- /dev/null +++ b/valida_shell.py @@ -0,0 +1,252 @@ +# Validazione del solutore shell su casi con soluzione nota. +# +# Non è una suite di test automatica: stampa i risultati e i confronti perché +# molti sono verifiche di convergenza, dove conta l'andamento dell'errore più +# di una soglia di pass/fail. Va eseguita a mano dopo ogni modifica alla +# fisica o all'assemblaggio: +# +# python valida_shell.py + +import math + +import numpy as np + +import assemblaggio as asm +import solutore as sol +import sorgente as src +from materials import MATERIALI +from mesh import campo_su_griglia, costruisci_mesh_cilindrica + + +MATERIALE = MATERIALI["banda_stagnata"] +K_TERM = MATERIALE["conducibilita_termica_W_mK"] +RHO = MATERIALE["densita_kg_m3"] +CP = MATERIALE["calore_specifico_J_kgK"] +ALPHA = K_TERM / (RHO * CP) + +LUNGHEZZA_M = 0.100 +RAGGIO_M = 0.035 +SPESSORE_M = 0.18e-3 +DT_S = 0.01 + + +def _esito(condizione: bool) -> str: + return "OK" if condizione else "FALLITO" + + +def _campo_nullo(mesh: dict) -> np.ndarray: + return np.zeros(mesh["n_nodi"]) + + +def verifica_matrici(mesh: dict) -> None: + print("=== 1. proprietà delle matrici assemblate ===") + C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M) + K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M) + + attesa = RHO * CP * SPESSORE_M * mesh["area_totale_m2"] + errore = abs(C.sum() / attesa - 1.0) + print(f" capacità totale {C.sum():.6f} vs rho*cp*s*A {attesa:.6f}" + f" err={errore:.1e} {_esito(errore < 1e-12)}") + + # A temperatura costante il flusso deve essere nullo, quindi ogni riga + # della matrice di conduzione somma a zero. + somma_righe = np.abs(np.asarray(K.sum(axis=1))).max() + print(f" max |somma righe K_cond| = {somma_righe:.1e}" + f" {_esito(somma_righe < 1e-9 * K_TERM)}") + + asimmetria = abs(K - K.T).max() + print(f" simmetria di K_cond: scarto {asimmetria:.1e}" + f" {_esito(asimmetria < 1e-12 * abs(K).max())}") + + autovalori = np.linalg.eigvalsh(K.toarray()[:400, :400]) + print(f" blocco campione semidefinito positivo:" + f" autovalore minimo {autovalori.min():.1e}" + f" {_esito(autovalori.min() > -1e-9)}") + + +def verifica_raffreddamento(mesh: dict) -> None: + # Campo uniforme senza sorgente: deve restare uniforme e decadere verso la + # temperatura ambiente con la costante di tempo rho*cp*s/(h_est+h_int). + # I bordi assiali sono esclusi, altrimenti il decadimento non è uniforme. + print("\n=== 2. raffreddamento convettivo uniforme ===") + aria = { + "temperatura_ambiente_C": 20.0, + "h_esterno_W_m2K": 12.0, + "h_interno_W_m2K": 8.0, + "h_bordi_W_m2K": 0.0, + } + C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M) + K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M) + K_conv, f_ambiente = asm.assembla_convezione(mesh, aria, SPESSORE_M) + solutore = sol.costruisci_solutore(C, K + K_conv, DT_S) + + sovratemperatura = 200.0 + T = np.full(mesh["n_nodi"], aria["temperatura_ambiente_C"] + sovratemperatura) + n_passi = 3000 + for _ in range(n_passi): + T = sol.passo_implicito(solutore, T, f_ambiente, _campo_nullo(mesh)) + + t = n_passi * DT_S + tau = RHO * CP * SPESSORE_M / (aria["h_esterno_W_m2K"] + aria["h_interno_W_m2K"]) + attesa = aria["temperatura_ambiente_C"] + sovratemperatura * math.exp(-t / tau) + errore = abs(T.mean() - attesa) + # Eulero implicito è del primo ordine: l'errore atteso è dell'ordine di + # dt/(2*tau) sulla sovratemperatura residua. + tolleranza = 3.0 * (attesa - aria["temperatura_ambiente_C"]) * DT_S / (2 * tau) + + print(f" tau = {tau:.3f} s, t = {t:.1f} s") + print(f" numerica {T.mean():.4f} °C vs analitica {attesa:.4f} °C" + f" err={errore:.4f} °C {_esito(errore < tolleranza)}") + dispersione = T.max() - T.min() + print(f" campo ancora uniforme: max-min = {dispersione:.1e} °C" + f" {_esito(dispersione < 1e-9)}") + + +def _sorgente_prova() -> dict: + return { + "x_inizio_m": 0.5, + "x_fine_m": 0.5, + "offset_y_percorso_m": 0.0, + "velocita_m_s": -1.0, + "numero_sorgenti": 1, + "distanza_sorgenti_m": 0.0, + "sigma_punto_m": 0.012, + "flusso_termico_picco_W_m2": 1e6, + "efficienza_riscaldamento": 1.0, + "zero_dopo_fine": True, + } + + +def verifica_sorgente(mesh: dict) -> None: + print("\n=== 3. vettore della sorgente e conservazione dell'energia ===") + sorgente = _sorgente_prova() + assemblatore = asm.prepara_assemblatore_sorgente(mesh) + fattore_arco = src.fattore_circonferenziale( + sorgente, RAGGIO_M, assemblatore["arco_gauss_m"] + ) + + # Con la gaussiana interamente dentro la fascetta la potenza assemblata + # deve valere l'integrale analitico 2*pi*sigma^2*q_max; con la sorgente + # centrata su un bordo assiale deve valerne esattamente la metà. + potenza_analitica = ( + sorgente["flusso_termico_picco_W_m2"] + * 2.0 * math.pi * sorgente["sigma_punto_m"] ** 2 + ) + # x_riferimento(t) = x_sensore + x_inizio + v*t, con x_sensore = 0.05. + for t, descrizione, atteso in [ + (0.50, "centrata a metà fascetta", 1.0), + (0.55, "centrata sul bordo x = 0", 0.5), + ]: + _, flusso_x = src.profilo_flusso_x_W_m2( + sorgente, 0.05, t, assemblatore["x_gauss_m"] + ) + f = asm.assembla_sorgente(assemblatore, flusso_x, fattore_arco) + rapporto = f.sum() / potenza_analitica + print(f" {descrizione:26s}: {f.sum():9.4f} W / {potenza_analitica:.4f} W" + f" = {rapporto:.4f} (atteso {atteso:.1f})" + f" {_esito(abs(rapporto - atteso) < 1e-3)}") + + # Senza convezione tutta l'energia iniettata deve finire nel campo. + aria = { + "temperatura_ambiente_C": 20.0, + "h_esterno_W_m2K": 0.0, + "h_interno_W_m2K": 0.0, + "h_bordi_W_m2K": 0.0, + } + C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M) + K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M) + _, f_ambiente = asm.assembla_convezione(mesh, aria, SPESSORE_M) + solutore = sol.costruisci_solutore(C, K, DT_S) + + T = np.full(mesh["n_nodi"], 20.0) + energia_iniziale = float((C @ T).sum()) + energia_iniettata = 0.0 + t = 0.0 + for _ in range(150): + _, flusso_x = src.profilo_flusso_x_W_m2( + sorgente, 0.05, t, assemblatore["x_gauss_m"] + ) + f = asm.assembla_sorgente(assemblatore, flusso_x, fattore_arco) + energia_iniettata += f.sum() * DT_S + T = sol.passo_implicito(solutore, T, f_ambiente, f) + t += DT_S + + incremento = float((C @ T).sum()) - energia_iniziale + errore = abs(incremento / energia_iniettata - 1.0) + print(f" iniettata {energia_iniettata:.4f} J, accumulata {incremento:.4f} J" + f" err={errore:.1e} {_esito(errore < 1e-10)}") + + +def verifica_modi_circonferenziali(mesh: dict) -> None: + # Un modo sin(n*theta) su un cilindro adiabatico decade come + # exp(-alpha * n^2 * t / R^2): confronto diretto con l'analitico. + print("\n=== 4. modi sinusoidali circonferenziali ===") + C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M) + K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M) + solutore = sol.costruisci_solutore(C, K, DT_S) + theta = np.tile(mesh["theta_nodi_rad"], mesh["n_nodi_x"]) + + for n in (1, 2, 4): + T = 100.0 * np.sin(n * theta) + for _ in range(200): + T = sol.passo_implicito( + solutore, T, _campo_nullo(mesh), _campo_nullo(mesh) + ) + t = 200 * DT_S + attesa = 100.0 * math.exp(-ALPHA * n**2 / RAGGIO_M**2 * t) + print(f" n={n}: numerica {np.abs(T).max():8.4f} vs analitica {attesa:8.4f}" + f" err={abs(np.abs(T).max() / attesa - 1):.3%}") + + +def verifica_convergenza() -> None: + print("\n=== 5. convergenza al raffinamento della mesh (modo n=4) ===") + for n_theta in (30, 60, 110, 220): + mesh = costruisci_mesh_cilindrica(LUNGHEZZA_M, RAGGIO_M, 20, n_theta) + C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M) + K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M) + solutore = sol.costruisci_solutore(C, K, DT_S) + theta = np.tile(mesh["theta_nodi_rad"], mesh["n_nodi_x"]) + T = 100.0 * np.sin(4 * theta) + for _ in range(200): + T = sol.passo_implicito( + solutore, T, _campo_nullo(mesh), _campo_nullo(mesh) + ) + attesa = 100.0 * math.exp(-ALPHA * 16 / RAGGIO_M**2 * 200 * DT_S) + print(f" n_elementi_theta={n_theta:4d} (arco {mesh['lato_arco_m']*1e3:5.2f} mm):" + f" err={abs(np.abs(T).max() / attesa - 1):.4%}") + + +def verifica_periodicita(mesh: dict) -> None: + # Un impulso su theta = 0 deve diffondere in modo identico nei due versi: + # è la prova che la connettività riavvolge davvero la circonferenza. + print("\n=== 6. periodicità circonferenziale ===") + C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M) + K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M) + solutore = sol.costruisci_solutore(C, K, DT_S) + + T = np.full(mesh["n_nodi"], 20.0) + T[50 * mesh["n_elementi_theta"]] = 500.0 + for _ in range(50): + T = sol.passo_implicito(solutore, T, _campo_nullo(mesh), _campo_nullo(mesh)) + + griglia = campo_su_griglia(mesh, T) + verso_positivo = griglia[50, 1:6] + verso_negativo = griglia[50, -1:-6:-1] + scarto = np.abs(verso_positivo - verso_negativo).max() + print(f" verso theta+ : {np.array2string(verso_positivo, precision=4)}") + print(f" verso theta- : {np.array2string(verso_negativo, precision=4)}") + print(f" scarto massimo = {scarto:.1e} {_esito(scarto < 1e-9)}") + + +def main() -> None: + mesh = costruisci_mesh_cilindrica(LUNGHEZZA_M, RAGGIO_M, 100, 110) + verifica_matrici(mesh) + verifica_raffreddamento(mesh) + verifica_sorgente(mesh) + verifica_modi_circonferenziali(mesh) + verifica_convergenza() + verifica_periodicita(mesh) + + +if __name__ == "__main__": + main()