davide/simulatore-induttori
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Trasforma il modello termico da 1D a 2D sulla sezione della fascetta

Browse Source 
La lastra diventa una fascetta cilindrica (diametro, lunghezza, spessore)
di cui si simula la sezione rettangolare lunghezza x spessore, con origine
nel vertice in alto a sinistra: x = lunghezza, z = spessore, y = coordinata
circonferenziale non risolta (offset collassato in attenuazione gaussiana).

- Le sorgenti si muovono in direzione -x e producono un profilo di flusso
  q(x) su tutto il lato esterno; corretto il posizionamento del gruppo per
  velocita' negative
- Sensore a infrarossi con coordinata x e distanza dalla parete interna,
  misura la superficie interna senza contatto
- Volumi finiti 2D con Eulero implicito: matrice sparsa fattorizzata LU una
  volta per run (scipy.sparse.linalg.splu), a ogni passo solo i sistemi
  triangolari
- Convezione sui lati esterno/interno, bordi adiabatici in x
- Schema CSV invariato; metadata con geometria e griglia 2D

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This commit is contained in:
Davide Grilli
2026-07-05 19:20:47 +02:00
parent b40ee49528
commit bb492c37a6
6 changed files with 189 additions and 147 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simulate.py
+126 -107
View File
@@ -6,8 +6,10 @@ from copy import deepcopy
from pathlib import Path
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import splu
from config import ARIA, PIASTRA, RANDOMIZZAZIONE, SENSORE, SIMULAZIONE, SORGENTE
from config import ARIA, FASCETTA, RANDOMIZZAZIONE, SENSORE, SIMULAZIONE, SORGENTE
from materials import MATERIALI
@@ -34,31 +36,28 @@ def _spread_sorgenti_m(sorgente: dict) -> float:
return (numero_sorgenti - 1) * distanza
def _x_riferimento_iniziale_m(sorgente: dict) -> float:
# Posizione a t=0 della sorgente di indice 0 (quella più arretrata nel
# verso di marcia). x_inizio_m è la distanza dal sensore della sorgente
# più avanzata (che quindi lo raggiunge per prima).
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_sensore = sorgente["x_sensore_m"]
x_inizio = sorgente["x_inizio_m"]
if sorgente["velocita_m_s"] >= 0:
return (x_sensore - x_inizio) - spread
return (x_sensore + x_inizio) + spread
return (x_sensore_m - x_inizio) - spread
return x_sensore_m + x_inizio
def _x_riferimento_finale_m(sorgente: dict) -> float:
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 (che quindi lo
# supera per ultima).
x_sensore = sorgente["x_sensore_m"]
# 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 + x_fine
return x_sensore - x_fine
def x_sorgente_al_tempo(sorgente: dict, t_s: float) -> float:
return _x_riferimento_iniziale_m(sorgente) + sorgente["velocita_m_s"] * t_s
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:
@@ -67,18 +66,21 @@ def _intervallo_attivo(inizio: float, fine: float, v: float, x_m: float) -> bool
return fine <= x_m <= inizio
def flusso_termico_incidente_W_m2(sorgente: dict, t_s: float) -> tuple[float, float]:
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
# flusso_termico_efficace_W_m2 (somma dei contributi di tutte le sorgenti).
# il profilo di flusso termico efficace q(x) [W/m²] sul lato esterno,
# somma dei contributi di tutte le sorgenti attive.
#
# Il movimento è rappresentato con un'impronta gaussiana centrata su
# ciascuna sorgente in moto. Più sorgenti equidistanti si muovono insieme
# come un gruppo rigido: condividono velocità, sigma e flusso di picco, e
# sono sfalsate lungo x di un multiplo di "distanza_sorgenti_m". Il
# modello 1D vede solo il flusso lungo la linea che passa per il sensore
# fisso, sommando il contributo di tutte le sorgenti attive.
x_rif_iniziale = _x_riferimento_iniziale_m(sorgente)
x_rif_finale = _x_riferimento_finale_m(sorgente)
# 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)
@@ -88,7 +90,10 @@ def flusso_termico_incidente_W_m2(sorgente: dict, t_s: float) -> tuple[float, fl
dy = sorgente["offset_y_percorso_m"]
sigma = sorgente["sigma_punto_m"]
q_totale = 0.0
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
@@ -98,74 +103,76 @@ def flusso_termico_incidente_W_m2(sorgente: dict, t_s: float) -> tuple[float, fl
if not _intervallo_attivo(inizio_i, fine_i, v, x_i):
continue
dx = x_i - sorgente["x_sensore_m"]
gaussiana = math.exp(-0.5 * (dx * dx + dy * dy) / (sigma * sigma))
q_totale += (
sorgente["flusso_termico_picco_W_m2"]
* sorgente["efficienza_riscaldamento"]
* gaussiana
)
dx = x_centri_m - x_i
q_x += q_picco * attenuazione_y * np.exp(-0.5 * (dx * dx) / (sigma * sigma))
return x_riferimento, q_totale
return x_riferimento, q_x
def riscaldamento_volumetrico_W_m3(
q_superficie_W_m2: float,
def profilo_deposizione_z_1_m(
z_centri_m: np.ndarray,
spessore_m: float,
skin_depth_m: float,
) -> np.ndarray:
# Converte il flusso superficiale equivalente in riscaldamento volumetrico q_vol(z).
# Profilo di deposizione volumetrica del flusso superficiale [1/m]:
#
# q_vol(z) = A * exp(-z / delta)
# p(z) = exp(-z / delta) / (delta * (1 - exp(-spessore / delta)))
#
# A è scelto in modo che integrale_0^L q_vol(z) dz = q_superficie_W_m2.
if q_superficie_W_m2 <= 0.0:
return np.zeros_like(z_centri_m)
# 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 (q_superficie_W_m2 / normalizzazione) * np.exp(-z_centri_m / delta)
return np.exp(-z_centri_m / delta) / normalizzazione
def costruisci_matrice_implicita(
n: int,
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,
materiale: dict,
h_caldo_W_m2K: float,
h_freddo_W_m2K: float,
) -> np.ndarray:
# Costruisce la matrice A per Eulero implicito:
h_esterno_W_m2K: float,
h_interno_W_m2K: float,
):
# Costruisce e fattorizza (LU sparsa) la matrice A per Eulero implicito 2D:
# A * T_next = rhs
#
# Le celle di bordo includono la convezione verso l'ambiente.
# Le incognite sono i centri cella T[i, j] con i lungo x e j lungo z,
# appiattiti in ordine C (indice = i * n_z + j). I bordi x = 0 e
# x = lunghezza sono adiabatici; i bordi z = 0 (esterno) e z = spessore
# (interno) includono la convezione verso l'ambiente.
k = materiale["conducibilita_termica_W_mK"]
rho = materiale["densita_kg_m3"]
cp = materiale["calore_specifico_J_kgK"]
alpha = k / (rho * cp)
r = alpha * dt_s / (dz_m * dz_m)
b_caldo = h_caldo_W_m2K * dt_s / (rho * cp * dz_m)
b_freddo = h_freddo_W_m2K * dt_s / (rho * cp * dz_m)
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)
A = np.zeros((n, n), dtype=float)
n = n_x * n_z
convezione = np.zeros(n)
convezione[0::n_z] = b_esterno
convezione[n_z - 1::n_z] = b_interno
# Bordo caldo, z = 0.
A[0, 0] = 1.0 + r + b_caldo
A[0, 1] = -r
# Celle interne.
for i in range(1, n - 1):
A[i, i - 1] = -r
A[i, i] = 1.0 + 2.0 * r
A[i, i + 1] = -r
# Bordo freddo, z = spessore.
A[n - 1, n - 2] = -r
A[n - 1, n - 1] = 1.0 + r + b_freddo
return A
A = (
sp.identity(n)
+ r_x * sp.kron(_laplaciano_1d(n_x), sp.identity(n_z))
+ r_z * sp.kron(sp.identity(n_x), _laplaciano_1d(n_z))
+ sp.diags(convezione)
)
return splu(sp.csc_matrix(A))
def quantizza(valore: float, passo: float) -> float:
@@ -175,7 +182,7 @@ def quantizza(valore: float, passo: float) -> float:
def configurazione_randomizzata(indice_run: int, rng: random.Random) -> dict:
piastra = deepcopy(PIASTRA)
fascetta = deepcopy(FASCETTA)
aria = deepcopy(ARIA)
sorgente = deepcopy(SORGENTE)
sensore = deepcopy(SENSORE)
@@ -214,7 +221,7 @@ def configurazione_randomizzata(indice_run: int, rng: random.Random) -> dict:
return {
"id_run": f"run_{indice_run:04d}",
"piastra": piastra,
"fascetta": fascetta,
"aria": aria,
"sorgente": sorgente,
"sensore": sensore,
@@ -222,18 +229,25 @@ def configurazione_randomizzata(indice_run: int, rng: random.Random) -> dict:
def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
piastra = cfg_run["piastra"]
fascetta = cfg_run["fascetta"]
aria = cfg_run["aria"]
sorgente = cfg_run["sorgente"]
sensore = cfg_run["sensore"]
nome_materiale = piastra["materiale"]
nome_materiale = fascetta["materiale"]
materiale = MATERIALI[nome_materiale]
spessore = piastra["spessore_mm"] / 1000.0
n = piastra["n_nodi"]
dz = spessore / n
z_centri = (np.arange(n) + 0.5) * dz
lunghezza = fascetta["lunghezza_mm"] / 1000.0
spessore = fascetta["spessore_mm"] / 1000.0
n_x = fascetta["n_nodi_x"]
n_z = fascetta["n_nodi_z"]
dx = lunghezza / n_x
dz = spessore / n_z
x_centri = (np.arange(n_x) + 0.5) * dx
z_centri = (np.arange(n_z) + 0.5) * dz
x_sensore = sensore["x_mm"] / 1000.0
i_sensore = min(n_x - 1, max(0, int(x_sensore / dx)))
dt = SIMULAZIONE["dt_interno_s"]
durata = SIMULAZIONE["durata_s"]
@@ -244,26 +258,30 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
else:
skin_depth = float(sorgente["skin_depth_fissa_m"])
A = costruisci_matrice_implicita(
n=n,
solutore = costruisci_solutore_implicito_2d(
n_x=n_x,
n_z=n_z,
dt_s=dt,
dx_m=dx,
dz_m=dz,
materiale=materiale,
h_caldo_W_m2K=aria["h_caldo_W_m2K"],
h_freddo_W_m2K=aria["h_freddo_W_m2K"],
h_esterno_W_m2K=aria["h_esterno_W_m2K"],
h_interno_W_m2K=aria["h_interno_W_m2K"],
)
rho = materiale["densita_kg_m3"]
cp = materiale["calore_specifico_J_kgK"]
b_caldo = aria["h_caldo_W_m2K"] * dt / (rho * cp * dz)
b_freddo = aria["h_freddo_W_m2K"] * dt / (rho * cp * dz)
b_esterno = aria["h_esterno_W_m2K"] * dt / (rho * cp * dz)
b_interno = aria["h_interno_W_m2K"] * dt / (rho * cp * dz)
T = np.full(n, piastra["temperatura_iniziale_C"], dtype=float)
T_sensore = T[-1]
profilo_z = profilo_deposizione_z_1_m(z_centri, spessore, skin_depth)
T = np.full((n_x, n_z), fascetta["temperatura_iniziale_C"], dtype=float)
T_sensore = T[i_sensore, -1]
prossimo_campione_t = 0.0
T_vera_max = T[-1]
T_vera_max = T[i_sensore, -1]
T_misurata_max = T_sensore
output_csv.parent.mkdir(parents=True, exist_ok=True)
@@ -289,22 +307,19 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
t = 0.0
while t <= durata + 1e-12:
x_sorgente, q_superficie = flusso_termico_incidente_W_m2(sorgente, t)
q_vol = riscaldamento_volumetrico_W_m3(
q_superficie_W_m2=q_superficie,
z_centri_m=z_centri,
spessore_m=spessore,
skin_depth_m=skin_depth,
x_sorgente, q_x = profilo_flusso_incidente_W_m2(
sorgente, x_sensore, t, x_centri
)
rhs = T + dt * q_vol / (rho * cp)
rhs[0] += b_caldo * aria["temperatura_ambiente_C"]
rhs[-1] += b_freddo * aria["temperatura_ambiente_C"]
rhs = T + (dt / (rho * cp)) * q_x[:, None] * profilo_z[None, :]
rhs[:, 0] += b_esterno * aria["temperatura_ambiente_C"]
rhs[:, -1] += b_interno * aria["temperatura_ambiente_C"]
T = np.linalg.solve(A, rhs)
T = solutore.solve(rhs.ravel()).reshape(n_x, n_z)
# Temperatura vera sul lato freddo, dove si trova il sensore fisso.
T_vera_lato_sensore = T[-1]
# Temperatura vera della superficie interna nel punto osservato
# dal sensore infrarosso.
T_vera_lato_sensore = T[i_sensore, -1]
# Inerzia del sensore del primo ordine.
tau_sensore = max(sensore["costante_tempo_s"], 1e-9)
@@ -323,11 +338,11 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
f"{t:.6f}",
f"{x_sorgente:.9f}",
f"{sorgente['offset_y_percorso_m']:.9f}",
f"{q_superficie:.6f}",
f"{q_x[i_sensore]:.6f}",
f"{skin_depth:.9e}",
f"{T_vera_lato_sensore:.6f}",
f"{misurata:.6f}",
f"{T[0]:.6f}",
f"{T[i_sensore, 0]:.6f}",
f"{aria['temperatura_ambiente_C']:.6f}",
f"{sorgente['velocita_m_s']:.9f}",
f"{sorgente['sigma_punto_m']:.9f}",
@@ -342,17 +357,21 @@ 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,
"lunghezza_m": lunghezza,
"spessore_m": spessore,
"n_nodi": n,
"n_nodi_x": n_x,
"n_nodi_z": n_z,
"durata_s": durata,
"frequenza_campionamento_hz": SIMULAZIONE["frequenza_campionamento_hz"],
"dt_interno_s": dt,
"temperatura_ambiente_C": aria["temperatura_ambiente_C"],
"h_caldo_W_m2K": aria["h_caldo_W_m2K"],
"h_freddo_W_m2K": aria["h_freddo_W_m2K"],
"h_esterno_W_m2K": aria["h_esterno_W_m2K"],
"h_interno_W_m2K": aria["h_interno_W_m2K"],
"x_inizio_m": sorgente["x_inizio_m"],
"x_fine_m": sorgente["x_fine_m"],
"x_sensore_m": sorgente["x_sensore_m"],
"x_sensore_m": x_sensore,
"distanza_sensore_parete_m": sensore["distanza_parete_mm"] / 1000.0,
"offset_y_percorso_m": sorgente["offset_y_percorso_m"],
"velocita_m_s": sorgente["velocita_m_s"],
"numero_sorgenti": sorgente.get("numero_sorgenti", 1),