357 lines
11 KiB
Python
357 lines
11 KiB
Python
|
import csv
|
||
|
|
import math
|
||
|
|
import random
|
||
|
|
import shutil
|
||
|
|
from copy import deepcopy
|
||
|
|
from pathlib import Path
|
||
|
|
|
||
|
|
import numpy as np
|
||
|
|
|
||
|
|
from config import ARIA, PIASTRA, RANDOMIZZAZIONE, SENSORE, SIMULAZIONE, SORGENTE
|
||
|
|
from materials import MATERIALI
|
||
|
|
|
||
|
|
|
||
|
|
MU0 = 4.0 * math.pi * 1e-7
|
||
|
|
|
||
|
|
|
||
|
|
def calcola_skin_depth_m(materiale: dict, frequenza_hz: float) -> float:
|
||
|
|
# Skin depth elettromagnetica approssimata:
|
||
|
|
# delta = sqrt(2 * rho_e / (omega * mu))
|
||
|
|
#
|
||
|
|
# Semplificata. Per acciai ferromagnetici il comportamento reale
|
||
|
|
# è fortemente non lineare con temperatura e campo magnetico.
|
||
|
|
rho_e = materiale["resistivita_elettrica_ohm_m"]
|
||
|
|
mu_r = materiale["permeabilita_relativa"]
|
||
|
|
omega = 2.0 * math.pi * frequenza_hz
|
||
|
|
mu = MU0 * mu_r
|
||
|
|
return math.sqrt(2.0 * rho_e / (omega * mu))
|
||
|
|
|
||
|
|
|
||
|
|
def x_sorgente_al_tempo(sorgente: dict, t_s: float) -> float:
|
||
|
|
return sorgente["x_inizio_m"] + sorgente["velocita_m_s"] * t_s
|
||
|
|
|
||
|
|
|
||
|
|
def sorgente_attiva(sorgente: dict, x_m: float) -> bool:
|
||
|
|
if not sorgente.get("zero_dopo_fine", True):
|
||
|
|
return True
|
||
|
|
|
||
|
|
inizio = sorgente["x_inizio_m"]
|
||
|
|
fine = sorgente["x_fine_m"]
|
||
|
|
v = sorgente["velocita_m_s"]
|
||
|
|
|
||
|
|
if v >= 0:
|
||
|
|
return inizio <= x_m <= fine
|
||
|
|
return fine <= x_m <= inizio
|
||
|
|
|
||
|
|
|
||
|
|
def flusso_termico_incidente_W_m2(sorgente: dict, t_s: float) -> tuple[float, float]:
|
||
|
|
# Restituisce x_sorgente_m e flusso_termico_efficace_W_m2.
|
||
|
|
#
|
||
|
|
# Il movimento è rappresentato con un'impronta gaussiana centrata sulla
|
||
|
|
# sorgente in moto. Il modello 1D vede solo il flusso lungo la linea
|
||
|
|
# che passa per il sensore fisso.
|
||
|
|
x = x_sorgente_al_tempo(sorgente, t_s)
|
||
|
|
|
||
|
|
if not sorgente_attiva(sorgente, x):
|
||
|
|
return x, 0.0
|
||
|
|
|
||
|
|
dx = x - sorgente["x_sensore_m"]
|
||
|
|
dy = sorgente["offset_y_percorso_m"]
|
||
|
|
sigma = sorgente["sigma_punto_m"]
|
||
|
|
|
||
|
|
gaussiana = math.exp(-0.5 * (dx * dx + dy * dy) / (sigma * sigma))
|
||
|
|
|
||
|
|
q = (
|
||
|
|
sorgente["flusso_termico_picco_W_m2"]
|
||
|
|
* sorgente["efficienza_riscaldamento"]
|
||
|
|
* gaussiana
|
||
|
|
)
|
||
|
|
return x, q
|
||
|
|
|
||
|
|
|
||
|
|
def riscaldamento_volumetrico_W_m3(
|
||
|
|
q_superficie_W_m2: float,
|
||
|
|
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).
|
||
|
|
#
|
||
|
|
# q_vol(z) = A * exp(-z / 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)
|
||
|
|
|
||
|
|
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)
|
||
|
|
|
||
|
|
|
||
|
|
def costruisci_matrice_implicita(
|
||
|
|
n: int,
|
||
|
|
dt_s: 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:
|
||
|
|
# A * T_next = rhs
|
||
|
|
#
|
||
|
|
# Le celle di bordo 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)
|
||
|
|
|
||
|
|
A = np.zeros((n, n), dtype=float)
|
||
|
|
|
||
|
|
# 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
|
||
|
|
|
||
|
|
|
||
|
|
def quantizza(valore: float, passo: float) -> float:
|
||
|
|
if passo <= 0.0:
|
||
|
|
return valore
|
||
|
|
return round(valore / passo) * passo
|
||
|
|
|
||
|
|
|
||
|
|
def configurazione_randomizzata(indice_run: int, rng: random.Random) -> dict:
|
||
|
|
piastra = deepcopy(PIASTRA)
|
||
|
|
aria = deepcopy(ARIA)
|
||
|
|
sorgente = deepcopy(SORGENTE)
|
||
|
|
sensore = deepcopy(SENSORE)
|
||
|
|
|
||
|
|
if RANDOMIZZAZIONE.get("abilitata", False):
|
||
|
|
def perturba_rel(valore: float, std_rel: float, fattore_min: float = 0.1) -> float:
|
||
|
|
fattore = rng.gauss(1.0, std_rel)
|
||
|
|
fattore = max(fattore_min, fattore)
|
||
|
|
return valore * fattore
|
||
|
|
|
||
|
|
sorgente["velocita_m_s"] = perturba_rel(
|
||
|
|
sorgente["velocita_m_s"],
|
||
|
|
RANDOMIZZAZIONE["velocita_std_rel"],
|
||
|
|
)
|
||
|
|
sorgente["flusso_termico_picco_W_m2"] = perturba_rel(
|
||
|
|
sorgente["flusso_termico_picco_W_m2"],
|
||
|
|
RANDOMIZZAZIONE["flusso_picco_std_rel"],
|
||
|
|
)
|
||
|
|
sorgente["sigma_punto_m"] = perturba_rel(
|
||
|
|
sorgente["sigma_punto_m"],
|
||
|
|
RANDOMIZZAZIONE["sigma_punto_std_rel"],
|
||
|
|
)
|
||
|
|
sorgente["offset_y_percorso_m"] = rng.uniform(
|
||
|
|
-RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
|
||
|
|
RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
|
||
|
|
)
|
||
|
|
aria["temperatura_ambiente_C"] += rng.gauss(
|
||
|
|
0.0,
|
||
|
|
RANDOMIZZAZIONE["temperatura_ambiente_std_C"],
|
||
|
|
)
|
||
|
|
sensore["rumore_std_C"] = perturba_rel(
|
||
|
|
sensore["rumore_std_C"],
|
||
|
|
RANDOMIZZAZIONE["rumore_sensore_std_rel"],
|
||
|
|
fattore_min=0.0,
|
||
|
|
)
|
||
|
|
|
||
|
|
return {
|
||
|
|
"id_run": f"run_{indice_run:04d}",
|
||
|
|
"piastra": piastra,
|
||
|
|
"aria": aria,
|
||
|
|
"sorgente": sorgente,
|
||
|
|
"sensore": sensore,
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
|
||
|
|
piastra = cfg_run["piastra"]
|
||
|
|
aria = cfg_run["aria"]
|
||
|
|
sorgente = cfg_run["sorgente"]
|
||
|
|
sensore = cfg_run["sensore"]
|
||
|
|
|
||
|
|
nome_materiale = piastra["materiale"]
|
||
|
|
materiale = MATERIALI[nome_materiale]
|
||
|
|
|
||
|
|
spessore = piastra["spessore_m"]
|
||
|
|
n = piastra["n_nodi"]
|
||
|
|
dz = spessore / n
|
||
|
|
z_centri = (np.arange(n) + 0.5) * dz
|
||
|
|
|
||
|
|
dt = SIMULAZIONE["dt_interno_s"]
|
||
|
|
durata = SIMULAZIONE["durata_s"]
|
||
|
|
periodo_campionamento = 1.0 / SIMULAZIONE["frequenza_campionamento_hz"]
|
||
|
|
|
||
|
|
if sorgente["skin_depth_fissa_m"] is None:
|
||
|
|
skin_depth = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
|
||
|
|
else:
|
||
|
|
skin_depth = float(sorgente["skin_depth_fissa_m"])
|
||
|
|
|
||
|
|
A = costruisci_matrice_implicita(
|
||
|
|
n=n,
|
||
|
|
dt_s=dt,
|
||
|
|
dz_m=dz,
|
||
|
|
materiale=materiale,
|
||
|
|
h_caldo_W_m2K=aria["h_caldo_W_m2K"],
|
||
|
|
h_freddo_W_m2K=aria["h_freddo_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)
|
||
|
|
|
||
|
|
T = np.full(n, piastra["temperatura_iniziale_C"], dtype=float)
|
||
|
|
T_sensore = T[-1]
|
||
|
|
|
||
|
|
prossimo_campione_t = 0.0
|
||
|
|
T_vera_max = T[-1]
|
||
|
|
T_misurata_max = T_sensore
|
||
|
|
|
||
|
|
output_csv.parent.mkdir(parents=True, exist_ok=True)
|
||
|
|
|
||
|
|
with output_csv.open("w", newline="") as f:
|
||
|
|
writer = csv.writer(f)
|
||
|
|
writer.writerow([
|
||
|
|
"id_run",
|
||
|
|
"tempo_s",
|
||
|
|
"x_sorgente_m",
|
||
|
|
"offset_y_sorgente_m",
|
||
|
|
"flusso_termico_sorgente_W_m2",
|
||
|
|
"skin_depth_m",
|
||
|
|
"T_vera_lato_sensore_C",
|
||
|
|
"T_misurata_sensore_C",
|
||
|
|
"T_lato_caldo_C",
|
||
|
|
"T_ambiente_C",
|
||
|
|
"velocita_m_s",
|
||
|
|
"sigma_punto_m",
|
||
|
|
"flusso_picco_W_m2",
|
||
|
|
"materiale",
|
||
|
|
])
|
||
|
|
|
||
|
|
t = 0.0
|
||
|
|
while t <= durata + 1e-12:
|
||
|
|
x_sorgente, q_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,
|
||
|
|
)
|
||
|
|
|
||
|
|
rhs = T + dt * q_vol / (rho * cp)
|
||
|
|
rhs[0] += b_caldo * aria["temperatura_ambiente_C"]
|
||
|
|
rhs[-1] += b_freddo * aria["temperatura_ambiente_C"]
|
||
|
|
|
||
|
|
T = np.linalg.solve(A, rhs)
|
||
|
|
|
||
|
|
# Temperatura vera sul lato freddo, dove si trova il sensore fisso.
|
||
|
|
T_vera_lato_sensore = T[-1]
|
||
|
|
|
||
|
|
# Inerzia del sensore del primo ordine.
|
||
|
|
tau_sensore = max(sensore["costante_tempo_s"], 1e-9)
|
||
|
|
T_sensore += (T_vera_lato_sensore - T_sensore) * dt / tau_sensore
|
||
|
|
|
||
|
|
# Campionamento CSV.
|
||
|
|
if t + 1e-12 >= prossimo_campione_t:
|
||
|
|
misurata = T_sensore + rng.gauss(0.0, sensore["rumore_std_C"])
|
||
|
|
misurata = quantizza(misurata, sensore["quantizzazione_C"])
|
||
|
|
|
||
|
|
T_vera_max = max(T_vera_max, T_vera_lato_sensore)
|
||
|
|
T_misurata_max = max(T_misurata_max, misurata)
|
||
|
|
|
||
|
|
writer.writerow([
|
||
|
|
cfg_run["id_run"],
|
||
|
|
f"{t:.6f}",
|
||
|
|
f"{x_sorgente:.9f}",
|
||
|
|
f"{sorgente['offset_y_percorso_m']:.9f}",
|
||
|
|
f"{q_superficie:.6f}",
|
||
|
|
f"{skin_depth:.9e}",
|
||
|
|
f"{T_vera_lato_sensore:.6f}",
|
||
|
|
f"{misurata:.6f}",
|
||
|
|
f"{T[0]:.6f}",
|
||
|
|
f"{aria['temperatura_ambiente_C']:.6f}",
|
||
|
|
f"{sorgente['velocita_m_s']:.9f}",
|
||
|
|
f"{sorgente['sigma_punto_m']:.9f}",
|
||
|
|
f"{sorgente['flusso_termico_picco_W_m2']:.6f}",
|
||
|
|
nome_materiale,
|
||
|
|
])
|
||
|
|
prossimo_campione_t += periodo_campionamento
|
||
|
|
|
||
|
|
t += dt
|
||
|
|
|
||
|
|
return {
|
||
|
|
"id_run": cfg_run["id_run"],
|
||
|
|
"file_csv": str(output_csv.name),
|
||
|
|
"materiale": nome_materiale,
|
||
|
|
"spessore_m": spessore,
|
||
|
|
"n_nodi": n,
|
||
|
|
"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"],
|
||
|
|
"x_inizio_m": sorgente["x_inizio_m"],
|
||
|
|
"x_fine_m": sorgente["x_fine_m"],
|
||
|
|
"x_sensore_m": sorgente["x_sensore_m"],
|
||
|
|
"offset_y_percorso_m": sorgente["offset_y_percorso_m"],
|
||
|
|
"velocita_m_s": sorgente["velocita_m_s"],
|
||
|
|
"sigma_punto_m": sorgente["sigma_punto_m"],
|
||
|
|
"flusso_termico_picco_W_m2": sorgente["flusso_termico_picco_W_m2"],
|
||
|
|
"efficienza_riscaldamento": sorgente["efficienza_riscaldamento"],
|
||
|
|
"frequenza_hz": sorgente["frequenza_hz"],
|
||
|
|
"skin_depth_m": skin_depth,
|
||
|
|
"costante_tempo_sensore_s": sensore["costante_tempo_s"],
|
||
|
|
"rumore_std_sensore_C": sensore["rumore_std_C"],
|
||
|
|
"quantizzazione_sensore_C": sensore["quantizzazione_C"],
|
||
|
|
"T_vera_max_lato_sensore_C": T_vera_max,
|
||
|
|
"T_misurata_max_sensore_C": T_misurata_max,
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
def main() -> None:
|
||
|
|
cartella_output = Path(SIMULAZIONE["cartella_output"])
|
||
|
|
if cartella_output.exists():
|
||
|
|
shutil.rmtree(cartella_output)
|
||
|
|
cartella_output.mkdir(parents=True, exist_ok=True)
|
||
|
|
|
||
|
|
rng = random.Random(SIMULAZIONE["seed"])
|
||
|
|
|
||
|
|
righe_metadata = []
|
||
|
|
for i in range(1, SIMULAZIONE["num_run"] + 1):
|
||
|
|
cfg_run = configurazione_randomizzata(i, rng)
|
||
|
|
percorso_csv = cartella_output / f"{cfg_run['id_run']}.csv"
|
||
|
|
righe_metadata.append(simula_singolo(cfg_run, percorso_csv, rng))
|
||
|
|
|
||
|
|
percorso_metadata = cartella_output / "metadata.csv"
|
||
|
|
with percorso_metadata.open("w", newline="") as f:
|
||
|
|
writer = csv.DictWriter(f, fieldnames=list(righe_metadata[0].keys()))
|
||
|
|
writer.writeheader()
|
||
|
|
writer.writerows(righe_metadata)
|
||
|
|
|
||
|
|
print(f"Generati {len(righe_metadata)} run in: {cartella_output.resolve()}")
|
||
|
|
print(f"Metadata: {percorso_metadata.resolve()}")
|
||
|
|
|
||
|
|
|
||
|
|
if __name__ == "__main__":
|
||
|
|
main()
|