davide
15510a06d1
Sostituisce la BC Neumann ghost-cell a x=0 con BC Robin su entrambi i bordi. Q(t) viene iniettato come termine sorgente puntuale al nodo più vicino a X_SRC, dopo le BCs per non essere sovrascritto. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
88 lines
2.3 KiB
Python
88 lines
2.3 KiB
Python
"""
|
|
FTCS (Forward-Time Centered-Space) explicit finite difference solver
|
|
for the 1D heat equation with a movable point heat source:
|
|
|
|
dT/dt = alpha * d²T/dx² + Q(t) * alpha / (k * dx) * delta(x - x_src)
|
|
|
|
Boundary conditions (both ends Robin / convective):
|
|
- x=0 (Robin): convective exchange with ambient
|
|
- x=L (Robin): convective exchange with ambient
|
|
|
|
Returns T_matrix of shape (NX, NT).
|
|
"""
|
|
|
|
import sys
|
|
import os
|
|
import numpy as np
|
|
|
|
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
|
|
import config
|
|
|
|
|
|
def solve():
|
|
"""Run the FTCS solver and return (T_matrix, x_vals, t_vals).
|
|
|
|
Returns
|
|
-------
|
|
T_matrix : np.ndarray, shape (NX, NT)
|
|
x_vals : np.ndarray, shape (NX,)
|
|
t_vals : np.ndarray, shape (NT,)
|
|
"""
|
|
alpha = config.ALPHA
|
|
k = config.K
|
|
L = config.L
|
|
T0 = config.T0
|
|
Q_val = config.Q_VAL
|
|
t_step = config.T_STEP
|
|
h = config.H_CONV
|
|
T_amb = config.T_AMB
|
|
T_end = config.T_END
|
|
NX = config.NX
|
|
NT = config.NT
|
|
x_src = config.X_SRC
|
|
|
|
x_vals = np.linspace(0.0, L, NX)
|
|
t_vals = np.linspace(0.0, T_end, NT)
|
|
dx = x_vals[1] - x_vals[0]
|
|
dt = t_vals[1] - t_vals[0]
|
|
|
|
r = alpha * dt / dx**2
|
|
if r > 0.5:
|
|
print(
|
|
f"[FDM WARNING] Stability condition violated: "
|
|
f"r={r:.4f} > 0.5 (dt={dt:.6g}, dx={dx:.6g}). "
|
|
"Solution may diverge."
|
|
)
|
|
|
|
i_src = int(np.argmin(np.abs(x_vals - x_src)))
|
|
|
|
robin_coeff = dx * h / k
|
|
|
|
T_matrix = np.zeros((NX, NT), dtype=np.float64)
|
|
T_matrix[:, 0] = T0
|
|
|
|
T_cur = np.full(NX, T0, dtype=np.float64)
|
|
|
|
for n in range(NT - 1):
|
|
t_now = t_vals[n]
|
|
|
|
T_new = T_cur.copy()
|
|
|
|
# Interior FTCS (nodi 1 .. NX-2)
|
|
T_new[1:-1] = T_cur[1:-1] + r * (T_cur[2:] - 2.0 * T_cur[1:-1] + T_cur[:-2])
|
|
|
|
# Robin BC a x=0
|
|
T_new[0] = (T_cur[1] + robin_coeff * T_amb) / (1.0 + robin_coeff)
|
|
|
|
# Robin BC a x=L
|
|
T_new[-1] = (T_cur[-2] + robin_coeff * T_amb) / (1.0 + robin_coeff)
|
|
|
|
# Sorgente puntuale al nodo i_src — applicata dopo le BCs per non essere sovrascritta
|
|
Q_now = Q_val if t_now >= t_step else 0.0
|
|
T_new[i_src] += Q_now * alpha * dt / (k * dx)
|
|
|
|
T_cur = T_new
|
|
T_matrix[:, n + 1] = T_cur
|
|
|
|
return T_matrix, x_vals, t_vals
|