From 15510a06d1b2985d0a8e69cd7ab921a0c4249891 Mon Sep 17 00:00:00 2001 From: Davide Grilli Date: Wed, 13 May 2026 22:02:22 +0200 Subject: [PATCH] generalizzazione sorgente FDM a posizione arbitraria x_src MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- fdm/solver.py | 88 +++++++++++++++++++++------------------------------ 1 file changed, 36 insertions(+), 52 deletions(-) diff --git a/fdm/solver.py b/fdm/solver.py index f703e17..d279e77 100644 --- a/fdm/solver.py +++ b/fdm/solver.py @@ -1,12 +1,12 @@ """ FTCS (Forward-Time Centered-Space) explicit finite difference solver -for the 1D heat equation: +for the 1D heat equation with a movable point heat source: - dT/dt = alpha * d²T/dx² + dT/dt = alpha * d²T/dx² + Q(t) * alpha / (k * dx) * delta(x - x_src) -Boundary conditions: - - x=0 (Neumann): heat flux step Q(t) applied via ghost cell - - x=L (Robin): convective boundary condition +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). """ @@ -15,7 +15,6 @@ import sys import os import numpy as np -# Allow importing config from the project root sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) import config @@ -26,76 +25,61 @@ def solve(): Returns ------- T_matrix : np.ndarray, shape (NX, NT) - Temperature at each spatial node (row) and time step (column). - T_matrix[i, n] = T(x_i, t_n). - x_vals : np.ndarray, shape (NX,) - Spatial node positions from 0 to L. - t_vals : np.ndarray, shape (NT,) - Time values from 0 to T_END. + x_vals : np.ndarray, shape (NX,) + t_vals : np.ndarray, shape (NT,) """ - # ----------------------------------------------------------------- - # Parameters from config - # ----------------------------------------------------------------- - alpha = config.ALPHA - k = config.K - L = config.L - T0 = config.T0 - Q_val = config.Q_VAL + 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 + 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 - # ----------------------------------------------------------------- - # Grid - # ----------------------------------------------------------------- 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] - # ----------------------------------------------------------------- - # Stability check (CFL for explicit heat equation: r <= 0.5) - # ----------------------------------------------------------------- r = alpha * dt / dx**2 - if dt > dx**2 / (2.0 * alpha): + if r > 0.5: print( f"[FDM WARNING] Stability condition violated: " - f"dt={dt:.6g} > dx²/(2*alpha)={dx**2/(2.0*alpha):.6g} (r={r:.4f} > 0.5). " + f"r={r:.4f} > 0.5 (dt={dt:.6g}, dx={dx:.6g}). " "Solution may diverge." ) - # ----------------------------------------------------------------- - # Allocate output matrix and set initial condition - # ----------------------------------------------------------------- - T_matrix = np.zeros((NX, NT), dtype=np.float64) - T_matrix[:, 0] = T0 # uniform IC + 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 - # Working array for the current time level T_cur = np.full(NX, T0, dtype=np.float64) - # ----------------------------------------------------------------- - # Time integration - # ----------------------------------------------------------------- for n in range(NT - 1): t_now = t_vals[n] - # --- Neumann BC at x=0: ghost cell --- - # Q(t) = Q_val if t >= t_step else 0 - Q = Q_val if t_now >= t_step else 0.0 - T_ghost = T_cur[1] + 2.0 * dx * Q / k # ghost node T[-1] - - # --- Interior FTCS update (indices 1 .. NX-2) --- 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]) - # --- Apply Neumann BC at i=0 using ghost cell --- - T_new[0] = T_cur[0] + r * (T_cur[1] - 2.0 * T_cur[0] + T_ghost) + # Robin BC a x=0 + T_new[0] = (T_cur[1] + robin_coeff * T_amb) / (1.0 + robin_coeff) - # --- Apply Robin BC at x=L (explicit) --- - T_new[-1] = (T_cur[-2] + dx * h / k * T_amb) / (1.0 + dx * h / k) + # 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