generalizzazione sorgente FDM a posizione arbitraria x_src
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>
This commit is contained in:
+36
-52
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user