9  Forma de la ecuación

La ecuación diferencial de bio-calor de Pennes (1948) modela la transferencia de calor en tejidos biológicos, integrando efectos de conducción, perfusión sanguínea y metabolismo. Su forma general es:

\[ \rho c \frac{\partial T}{\partial t} = k_{\text{eff}} \frac{\partial^{2} T}{\partial x^{2}} - \rho_b c_b \omega_b (T - T_a) + \mathcal{Q}, \quad x \in \Omega, \, t \in [0, t_f] \tag{9.1}\]

Tabla 9.1: Tabla de nomenclaturas de la Ecuación 9.1.
Símbolo Descripción Unidades
\(T\) Temperatura del tejido °C
\(\rho\) Densidad del tejido \(\frac{kg}{m^3}\)
\(c\) Calor específico del tejido \(\frac{J}{kg °C}\)
\(k_{\text{eff}}\) Conductividad térmica \(\frac{W}{m °C}\)
\(\rho_b\) Densidad de la sangre \(\frac{kg}{m^3}\)
\(c_b\) Calor específico de la sangre \(\frac{J}{kg °C}\)
\(\omega_b\) Tasa de perfusión sanguínea \(1/s\)
\(T_a\) Temperatura arterial °C
\(\mathcal{Q} = q_m + q_p\) Fuente de calor \(\frac{W}{m^3}\)
\(q_m\) Metabolismo \(\frac{W}{m^3}\)
\(q_p\) Externa \(\frac{W}{m^3}\)

9.1 Versión reducida (adimensionalizada)

Mediante escalamiento: \[ T' = T - T_a \quad \theta = \dfrac{T'}{T_M - T_a} \quad X = \dfrac{x}{L_0} \quad \tau = \dfrac{t}{t_f} \tag{9.2}\]

Tabla 9.2: Tabla de nomenclatura de las relaciones para escalamiento.
Símbolo Descripción Unidades
\(L_0\) Longitud característica del dominio \(m\)
\(t_f\) Tiempo final de simulación \(s\)

la Ecuación 9.1 se convierte en:

\[ \partial_{\tau} \theta = a_1 \partial_{XX} \theta - a_2 W \theta + a_3 \] para una dimensión espacial; para el caso dos dimensional se tiene: \[ \partial_{\tau} \theta = a_1 \nabla^2 \theta - a_2 W \theta + a_3 \tag{9.3}\]

Parámetros adimensionales:
- \(a_1 = \frac{t_f}{\alpha L_0^2}\) (difusividad térmica \(\alpha = \frac{k_{\text{eff}}}{\rho c}\)).
- \(a_2 = \frac{t_f c_b}{\rho c}\).
- \(a_3 = \frac{t_f \mathcal{Q}}{\rho c (T_M - T_a)}\).
- \(W = \rho_b \omega_b\): Tasa volumétrica de perfusión (kg/m³·s).

Cabe decir que se modeló el caso más sencillo, que es asumiendo la fuente de calor \(\mathcal{Q} = 0\).

9.2 Condiciones de uso adecuadas

  1. Tejidos homogéneos: Aproximación válida para regiones con propiedades térmicas uniformes.
  2. Perfusión sanguínea constante: Supone flujo sanguíneo estable en el dominio.
  3. Aplicaciones clínicas: Hipertermia, crioterapia y modelado térmico en terapias oncológicas.

9.3 Solución analítica

Consideremos la ecuación diferencial parcial 9.3 sin el término \(a_3\), \[ \partial _\tau \theta(x,y,\tau) = a_1 \nabla^2(x,y,\tau) - a_2W\theta(x,y,\tau), \tag{9.4}\]

donde \(a_{1}\) y \(a_{2}\) son parámetros positivos adimensionales, \(W\) es una constante asociada al término de disipación. El dominio de estudio corresponde al cuadrado \([0,1]\times[0,1]\) en el espacio y al intervalo \([0,1]\) en el tiempo adimensional \(\tau\).

Las condiciones de frontera establecidas son mixtas:

  • En \(y=0\) es \(y=1\) se imponen condiciones de tipo Neumann, es decir, \[ ∂_y ​θ(x,0,τ)=0, \ \ ∂_y ​θ(x,1,τ)=0, \ \ τ≥0. \]

  • En \(x=0\) se prescribe una condición de tipo Dirichlet: \[ θ(0,y,τ)=0, \ \ τ≥0. \]

  • En \(x=1\) se fija una condición de tipo Neumann no homogénea, con dependencia lineal en el tiempo: \[ ∂_x ​θ(1,y,τ)=τ, \ \ τ≥0. \]

Además, se establece la condición inicial \[ θ(x,y,0)=0,\ \ (x,y)∈[0,1]×[0,1]. \]

9.3.1 Reducción del problema

Obsérvese que la ecuación es independiente de la variable \(y\), y que las condiciones de frontera en esa dirección son homogéneas (Neumann). Por tanto, la solución puede asumirse también independiente de \(y\), reduciendo el problema a una dimensión espacial. En consecuencia, se escribe \[ θ(x,y,τ) \equiv θ(x,τ). \]

La ecuación 9.4 se reduce a resolver \[ \partial _\tau \theta(x,\tau) = a_1 \dfrac{\partial^2 \theta}{\partial x^2}(x,\tau) - a_2W\theta(x,\tau), \tag{9.5}\]

con las condiciones de frontera \[ θ(0,τ)=0, \ \ ∂_x​θ(1,τ)=τ, \]

y la condición inicial \[ \theta(x,0) = 0. \]

9.3.2 Método de solución

Para resolver este problema, se aplica la técnica de separación en solución particular + solución homogénea. Se introduce la transformación

\[ θ(x,τ)=u(x,τ)+φ(x,τ), \]

donde \(\varphi(x,\tau)\) se escoge de manera que satisfaga la condición de frontera no homogénea en \(x=1\). Una elección natural es

\[ φ(x,τ)=xτ, \]

pues se cumple \[ ∂_x​φ(1,τ)=τ, \ \ φ(0,τ)=0. \]

De este modo, la función \(u(x,\tau)\) obedece condiciones homogéneas: \[ u(0,τ)=0, \ \ ∂_x​u(1,τ)=0, \]

y al sustituir en la ecuación 9.5, se obtiene para \(u(x,\tau)\): \[ ∂_τ​u = a_1\dfrac{​∂^2 u}{∂x^2}​ − a_2​Wu − a_2​Wxτ−x. \]

9.3.3 Solución mediante series

Se plantea una expansión de \(u(x,\tau)\) en términos de los autovalores y autofunciones del problema de Sturm–Liouville asociado: \[ u(x,τ)= \sum_{m=1}^{\infty}​a_m​(τ) sen(α_m​x), \]

donde los modos propios satisfacen \[ \alpha_m = \dfrac{(2m-1)\pi}{2}. \]

La evolución temporal de los coeficientes \(a_{m}(\tau)\) resulta \[ a_m(\tau) = c_m\dfrac{1-e^{-\sigma_m\tau}}{\sigma_m} + d_m\dfrac{\sigma_m \tau-1+e^{-\sigma_m\tau}}{\sigma_m^2}, \]

con \[ \sigma_m = a_1\alpha_m^2 + a_2W, \ \ c_m = -\dfrac{2(-1)^{m-1}}{\alpha_m^2}, \ \ d_m = -\dfrac{2a_2W(-1)^{m-1}}{\alpha_m^2}. \]

Finalmente la solución aproximada truncada a \(M\) modos se expresa como \[ u(x,τ) \approx x\tau + \sum_{m=1}^{M}​a_m​(τ)sen(α_m​x). \tag{9.6}\]

9.3.4 Solución truncada codificada

Basado la Ecuación 9.6 se codificó la solución analítica aproximada y se guardó en un dataframe para los tiempos de interés.

Listado 9.1: Guardado de los datos de la solución analítica.
Código 
import numpy as np
import pandas as pd

# ---------- Parámetros físicos ----------
p = 1050     # densidad
c = 3639     # calor específico
keff = 5     # conductividad efectiva
tf = 1800    # tiempo característico
L0 = 0.05    # longitud característica
cb = 3825    # perfusión
Q = 0        # fuente (no usada aquí)

# Parámetro auxiliar
alpha_phys = p * c / keff

# ---------- Coeficientes adimensionales ----------
a1 = tf / (alpha_phys * L0**2)
a2 = tf * cb / (p * c)

# ---------- Parámetros de la serie analítica ----------
W = 1.0
M = 60  # número de modos en el truncamiento

m = np.arange(1, M+1)
alpha_m = (2*m - 1) * np.pi / 2.0     # α_m
lambda_m = alpha_m**2
sigma_m = a1 * lambda_m + a2 * W      # σ_m

sign = (-1.0)**(m-1)
c_m = -2.0 * sign / (alpha_m**2)
d_m = -2.0 * a2 * W * sign / (alpha_m**2)

# ---------- Mallado en espacio y tiempos ----------
step = 0.04
# 0.00, 0.04, ..., 1.00
grid_vals = np.round(np.arange(0.0, 1.0 + 1e-12, step), 2)  
x = grid_vals.copy()
y = grid_vals.copy()
X, Y = np.meshgrid(x, y, indexing='xy')

times = [0.0, 0.25, 0.5, 0.75, 1.0]

# ---------- Función para calcular θ(x,τ) ----------
def theta_xt(x_vec, tau_val):
    """Devuelve theta(x, tau) en un vector de x,
    usando la serie truncada."""
    # coeficientes a_m(tau)
    small_mask = np.isclose(sigma_m, 0.0, atol=1e-12)
    a = np.empty_like(sigma_m)
    if np.any(small_mask):
        a[small_mask] = (c_m[small_mask]*tau_val 
        + 0.5*d_m[small_mask]*tau_val**2)
    if np.any(~small_mask):
        s = sigma_m[~small_mask]
        cm = c_m[~small_mask]
        dm = d_m[~small_mask]
        a[~small_mask] = (cm*(1.0 - np.exp(-s*tau_val))/s 
        + dm*(s*tau_val - 1.0 
        + np.exp(-s*tau_val))/s**2)
    # suma modal
    sin_ax = np.sin(np.outer(alpha_m, x_vec))   # (M, Nx)
    theta_series = a.dot(sin_ax)                # (Nx,)
    # añadir φ(x,τ) = x τ
    return theta_series + tau_val * x_vec

# ---------- Construcción del DataFrame ----------
results = []
for t_val in times:
    theta_x = theta_xt(x, t_val)         # (Nx,)
    Theta = np.tile(theta_x, (y.size, 1))  # shape (Ny, Nx)
    for xi, yi, thetai in zip(X.ravel(), Y.ravel(), Theta.ravel()):
        results.append([t_val, float(xi), float(yi), float(thetai)])

df = pd.DataFrame(results, columns=["time", "X", "Y", "Theta"])

# ---------- Guardar en CSV ----------
df.to_csv("data/sol_analitica.csv", index=False)