Análisis nodal por diferencias finitas en transferencia de calor: guía paso a paso

Interacción original de Ciencia Esencial.

En el bloque interactivo anterior observaste que un nodo interior, uno de frontera y uno de esquina no representan la misma porción del sólido ni están sometidos a los mismos mecanismos. La visualización mostró conducción, temperatura prescrita, convección, flujo de calor impuesto y aislamiento; sin embargo, la misma lógica también permite incorporar radiación, generación interna y geometrías en las que las áreas o las distancias entre nodos dejan de ser uniformes. En todos los casos, la formulación parte de una sola idea: aplicar conservación de la energía al volumen de control asociado al nodo.

En estado estacionario no existe acumulación de energía dentro del volumen nodal. Por tanto, la energía que entra, la que sale y la que se genera deben compensarse:

Forma compacta del balance de energía para un nodo en estado estacionario.

Al escribir de manera explícita los mecanismos que pueden actuar sobre el nodo P, el balance adopta una forma más amplia:

Balance nodal general: solo se conservan los términos correspondientes a los mecanismos que realmente actúan.

Qué representa cada símbolo

La ecuación (2) no debe memorizarse como una expresión rígida. Funciona como una estructura de trabajo: si una superficie está aislada, su flujo es cero; si no existe radiación, se elimina ese término; si el material no genera calor, desaparece q‴V_P. En una malla cartesiana uniforme, la distancia internodal Δ_j corresponde a Δx o Δy, según la dirección analizada. En una frontera irregular o una malla no uniforme, debe utilizarse la distancia real entre el nodo P y su vecino j.

Una vez establecido el balance general, el siguiente paso es definir cuánto material representa cada temperatura nodal. En una malla rectangular, los límites del volumen de control se sitúan a mitad de distancia entre nodos vecinos. Por esta razón, un nodo interior conserva la extensión completa de la celda, mientras que un nodo ubicado sobre una frontera representa solo la mitad en la dirección normal a esa superficie.

Para una profundidad unitaria perpendicular al plano, el volumen asociado al nodo P puede expresarse de forma compacta como:

Los factores geométricos f_x y f_y valen 1 o 1/2 según la posición del nodo; b es el espesor perpendicular al plano. Para un análisis por unidad de profundidad se adopta b = 1.

Estas fracciones no solo describen la geometría. También determinan la magnitud del término de generación q‴V_P y las áreas disponibles para conducción, convección o radiación. Por ello, una ecuación de frontera o de esquina no puede obtenerse copiando directamente la formulación de un nodo interior.

Las secciones anteriores ya establecieron el balance general y el volumen representado por cada nodo. Aquí aplicaremos esas ideas al caso más sencillo: un nodo interior de una placa homogénea, rodeado por cuatro vecinos y sin contacto directo con una frontera externa.

La orientación añadida a la malla interactiva permite identificar la notación que utilizaremos:

1.4.1 Nodo interior sin generación

Para conducción bidimensional en estado estacionario, con conductividad térmica constante y sin generación interna, la ecuación diferencial gobernante es:

Ecuación de Laplace para conducción bidimensional estacionaria.

Al aproximar las segundas derivadas mediante diferencias centrales en una malla rectangular, se obtiene:

Formulación válida cuando los espaciamientos horizontal y vertical pueden ser diferentes.

Si la malla es cuadrada, de modo que Δx = Δy = Δ, la ecuación se reduce a:

Ecuación algebraica del nodo interior sin generación.

Al despejar la temperatura del nodo analizado:

En ausencia de generación, la temperatura interior es el promedio de sus cuatro vecinos.

1.4.2 Nodo interior con generación uniforme

Si el material produce calor a una razón volumétrica uniforme q‴, la ecuación diferencial incluye un término adicional:

Ecuación de Poisson para un material con generación volumétrica uniforme.

Antes de simplificar el espesor de la placa, el balance energético del volumen interior se escribe con todas las áreas y el volumen reales:

Cada término contiene el espesor constante b. Al factorizarlo y dividir toda la igualdad entre b, el espesor se cancela sin modificar el balance.

Después de esta simplificación, para una malla rectangular la formulación discreta correspondiente es:

El término de generación modifica el balance aun cuando el nodo permanezca en el interior.

Para una malla cuadrada, la ecuación se simplifica a:

Ecuación algebraica del nodo interior con generación uniforme.

Al despejar la temperatura nodal:

La generación añade una contribución positiva a la temperatura del nodo cuando q‴ > 0.

2.1 Temperatura prescrita

Cuando la temperatura de una superficie es conocida, el nodo de frontera deja de ser una incógnita. Su valor se asigna directamente y se incorpora en las ecuaciones de los nodos vecinos.

Condición de Dirichlet: la temperatura superficial es un dato del problema.

En la visualización interactiva, la frontera superior se fijó en T_sup = 200 °C. Si esa temperatura corresponde al vecino norte de un nodo interior, la ecuación del nodo interior se escribe como:

La temperatura conocida se traslada al término independiente.

2.2 Flujo de calor prescrito

En esta condición se conoce la energía que atraviesa la superficie por unidad de área, mientras que la temperatura de la frontera permanece desconocida. Para el medio volumen asociado a un nodo sobre una frontera vertical, el flujo impuesto se suma al balance como q″A.

Nodo de frontera vertical con flujo superficial prescrito y generación uniforme.

El coeficiente 2 delante de T_E aparece porque la cara orientada hacia el interior posee el área completa Δb, mientras que las caras norte y sur del medio volumen tienen área Δb/2.

2.3 Superficie aislada y plano de simetría

Una superficie aislada es un caso particular de flujo prescrito en el que no existe transferencia de calor normal al contorno:

Condición adiabática o de simetría: el gradiente normal es nulo.

Al sustituir q″ = 0 en la ecuación (14), resulta:

La formulación también representa una línea de simetría térmica.

2.4 Frontera convectiva

Cuando la superficie intercambia calor con un fluido, el aporte convectivo depende de la diferencia entre la temperatura del fluido y la temperatura superficial:

El signo ya es coherente con la convención de energía positiva hacia el volumen nodal.

Para abreviar los coeficientes se define el número de Biot de malla:

Aquí la longitud característica es el espaciamiento nodal, no la dimensión global del sólido.

El balance del nodo convectivo queda:

La convección modifica tanto el coeficiente de T_P como el término independiente.

2.5 Frontera radiativa

En una superficie que intercambia radiación con unos alrededores grandes, la tasa neta que entra al volumen nodal puede escribirse como:

Las temperaturas radiativas deben expresarse en kelvin.

Al incorporar este término al balance del medio volumen:

La radiación introduce una dependencia no lineal con la temperatura nodal.

2.6 Convección, radiación y flujo combinados

Una superficie real puede recibir un flujo impuesto y, al mismo tiempo, intercambiar energía por convección y radiación. En ese caso no se requiere una nueva teoría: el balance se construye sumando únicamente los mecanismos presentes.

Balance combinado para una frontera vertical de una malla cuadrada.

2.7 Interfaz entre materiales

En una interfaz con contacto térmico perfecto, la temperatura es continua, pero los flujos de ambos lados deben calcularse con la conductividad térmica (k), el área normal a la conducción (A) y la distancia o espaciamiento internodal (Δ) correspondientes a cada material.

Cada conexión queda ponderada por su conductancia térmica kA/Δ.

Si las áreas y los espaciamientos son iguales y no existe generación, la expresión se reduce a:

Caso particular de contacto perfecto y generación nula.

2.8 Esquinas y combinación de condiciones

Una esquina representa un cuarto del volumen interior y posee dos superficies expuestas. Por ello, su ecuación se construye sumando los mecanismos de ambas caras y la conducción hacia sus dos vecinos interiores.

2.8.1 Esquina con doble convección

Para una esquina con convección en las superficies 1 y 2, vecinos interiores al oeste y al norte, y generación uniforme:

Forma normalizada directamente a partir del cuarto de volumen de control.

2.8.2 Otras combinaciones

2.9 Fronteras irregulares

Cuando el contorno corta la malla, las áreas y los volúmenes dejan de ser fracciones simples de Δ². La ecuación general del balance sigue siendo válida, pero debe utilizarse la geometría real del volumen de control.

A_j, Δ_j y V_P deben calcularse para la porción real del dominio.

Refinar la malla cerca del contorno mejora la representación geométrica, pero no corrige un balance mal construido. Antes de resolver, deben verificarse las áreas expuestas, las distancias internodales, el volumen recortado y los mecanismos que actúan sobre cada segmento de frontera.

3.1 Enunciado gráfico del problema

Se analizará una placa cuadrada de 0,20 m de lado y conductividad k = 20 W/(m·K). La frontera izquierda se mantiene a 100 °C y la superior a 200 °C. La frontera derecha intercambia calor por convección con un fluido a T_∞,r = 25 °C y h_r = 50 W/(m²·K), mientras que la frontera inferior lo hace con T_∞,b = 30 °C y h_b = 20 W/(m²·K). No existe generación interna y se emplea una malla de Δ = 0,10 m.

3.2 Qué ecuación corresponde a cada nodo

Antes de formular, conviene revisar qué condiciones de las secciones 1.4 y 2 actúan realmente en cada nodo desconocido. La tabla siguiente no resuelve el problema; su función es indicar qué modelo nodal debe buscarse en cada caso.

• T₁: se utiliza la ecuación del nodo interior sin generación, porque el nodo está completamente rodeado por sólido y solo recibe temperaturas conocidas en dos direcciones.
• T₂: se utiliza la ecuación de frontera convectiva. La orientación geométrica cambia respecto al esquema general, pero el modelo es el mismo; simplemente se rota hacia la frontera derecha.
• T₃: se vuelve a usar la ecuación de frontera convectiva, ahora rotada hacia la frontera inferior.
• T₄: se aplica la ecuación de esquina con doble convección, ya que el cuarto de volumen está expuesto en dos caras.

3.3 Planteamiento del sistema de ecuaciones

Como las fronteras derecha e inferior son convectivas, primero se calculan sus números de Biot de malla. Estos valores permiten sustituir la convección sin reescribir el coeficiente hΔ/k en cada paso.

El orden de las incógnitas será [T₁, T₂, T₃, T₄]. Ese orden debe mantenerse de forma consistente cuando se construya cada fila de la matriz.

3.3.1 Nodo interior T₁

Para T₁ se usa la ecuación del nodo interior sin generación. Al sustituir las fronteras conocidas, 100 °C al oeste y 200 °C al norte, la fila queda:

Esta ecuación incorpora directamente las dos temperaturas prescritas vecinas.

3.3.2 Frontera convectiva derecha T₂

Para T₂ se emplea la ecuación de frontera convectiva. El vecino interior es T₁, el nodo superior conocido vale 200 °C y el nodo inferior vecino es T₄:

La convección aparece en el coeficiente de T₂ y en el término independiente.

3.3.3 Frontera convectiva inferior T₃

Para T₃ se usa la misma ecuación convectiva, pero rotada hacia la frontera inferior. El vecino interior es T₁, la esquina izquierda conocida vale 100 °C y el vecino derecho es T₄:

3.3.4 Esquina con doble convección T₄

Para T₄ se aplica la ecuación de esquina con doble convección. El cuarto de volumen conduce hacia T₂ y T₃, y convectivamente intercambia calor con el ambiente derecho e inferior:

La última línea se obtiene multiplicando toda la ecuación por 2 para escribirla con coeficientes comparables a los de las demás filas.

3.3.5 Forma matricial y resolución

Al ordenar los coeficientes según el vector de incógnitas [T₁, T₂, T₃, T₄], el problema adopta la forma matricial siguiente:

En este punto, la solución puede obtenerse con una calculadora científica que resuelva sistemas lineales, en Excel, en GNU Octave o con software matemático y de programación.

• Calculadora científica: útil cuando el sistema es pequeño y el modelo ya está correctamente formulado.
• Excel o una hoja de cálculo: permite organizar la matriz y aplicar herramientas de álgebra lineal o fórmulas matriciales.
• GNU Octave, MATLAB o Python/SciPy: adecuados para sistemas medianos y para automatizar la construcción de la matriz.

3.4 Solución e interpretación física

3.5 Cuándo el sistema es lineal y cuándo deja de serlo

La solución obtenida arriba es válida porque el problema conduce a un sistema lineal. Sin embargo, no todas las formulaciones de diferencias finitas conservan esa propiedad. La tabla siguiente resume qué ocurre en la ecuación nodal y qué tipo de sistema resulta en cada caso.

4.1 Verificación de la solución

4.1.1 Coherencia física

La primera comprobación no requiere software: consiste en volver al enunciado y preguntarse si el campo térmico obtenido responde a las condiciones que lo originaron. Para el ejemplo de la sección 3, las temperaturas desconocidas se encuentran entre 110,19 °C y 137,43 °C, dentro del intervalo impuesto por las fronteras y los fluidos del entorno.

4.1.2 Residuos de las ecuaciones nodales

Después de resolver el sistema, cada temperatura se sustituye nuevamente en las ecuaciones nodales. La diferencia entre el lado izquierdo y el término conocido forma el vector de residuos:

Cada componente r_i corresponde al incumplimiento algebraico de una ecuación nodal.

Para resumir la comprobación puede utilizarse el mayor residuo absoluto:

Cuanto más próximo a cero sea este valor, mejor se satisface el sistema algebraico.

Con las temperaturas internas conservadas por el programa, los residuos del ejemplo son del orden de 10⁻¹⁴, es decir, próximos a la precisión numérica de la operación. Al sustituir las temperaturas publicadas con dos decimales aparecen residuos pequeños —aproximadamente 0,010; 0,005; 0,014 y −0,027— producidos por el redondeo.

4.1.3 Balance energético global

Esta comprobación vuelve a mirar la placa como un solo sistema. En estado estacionario, el calor que entra y el que se genera deben compensar el calor que sale. Para comprobarlo sin perderse entre nodos, conviene separar el cálculo en tres cantidades:

1. Calcula el aporte de cada frontera por separado, utilizando la temperatura y el área que realmente representa cada nodo.
2. Suma todos los aportes que entran, añade la generación interna y suma por otra parte todo lo que sale.
3. Sustituye esos tres totales en la ecuación (37). Un valor pequeño de ε_E indica que el balance global cierra de manera satisfactoria.
4. Si el error es grande, revisa primero signos, áreas de frontera y esquinas antes de culpar al método de solución.

La tolerancia aceptable debe definirse de acuerdo con la precisión buscada y la escala del problema.

Para sumar el calor de una cara discretizada, imagina que cada nodo controla un tramo de esa superficie. Un nodo situado entre dos extremos representa un tramo completo de longitud Δ; un nodo de esquina representa solo medio tramo, Δ/2, porque la otra mitad pertenece a la cara perpendicular. Así se evita contar dos veces la zona próxima a la esquina.

4.1.4 Independencia de malla

Los residuos pueden ser casi nulos incluso con una malla demasiado gruesa, porque únicamente miden el cumplimiento de las ecuaciones discretizadas. Para estimar el error de discretización se repite el problema con una malla más fina y se comparan magnitudes evaluadas en posiciones equivalentes.

La comparación debe realizarse en puntos que representen la misma posición física.

1. Resuelve el problema con un espaciamiento inicial Δ razonable y conserva las temperaturas y los flujos de interés.
2. Repite el cálculo con Δ/2, reconstruyendo todas las ecuaciones y no solo ampliando la matriz anterior.
3. Compara temperaturas, flujos superficiales y calor total en posiciones equivalentes. Los flujos suelen ser más sensibles al refinamiento que las temperaturas.
4. Si el cambio supera la tolerancia definida, repite con Δ/4 o aplica refinamiento local en regiones con gradientes elevados.
5. Acepta la malla cuando una nueva reducción del espaciamiento ya no modifique de forma significativa las magnitudes relevantes.

Reducir el espaciamiento disminuye normalmente el error de discretización, pero aumenta el número de operaciones y puede acumular más error de redondeo. Por esa razón, una comprobación adicional consiste en repetir el mismo cálculo con mayor precisión numérica y verificar que la solución no cambie de manera apreciable.

4.2 Implementación y comprobación con software

El software no decide qué volumen representa cada nodo ni qué condición debe aparecer en el balance. Esa parte corresponde al modelado desarrollado en las secciones anteriores. Su función comienza después: organizar los coeficientes, resolver el sistema sin perder precisión y automatizar las verificaciones.

4.2.1 Python con NumPy

El siguiente código reproduce el sistema del ejemplo, resuelve las cuatro temperaturas y calcula automáticamente los residuos. Se utiliza el nombre rhs para el vector del lado derecho y evitar confundirlo con el espesor b de la placa.

from __future__ import annotations

import numpy as np

A = np.array([
    [ 4.0, -1.0, -1.0,  0.0],
    [-2.0,  4.5,  0.0, -1.0],
    [-2.0,  0.0,  4.2, -1.0],
    [ 0.0, -2.0, -2.0,  4.7],
], dtype=float)

rhs = np.array([300.0, 212.5, 106.0, 18.5], dtype=float)

if A.ndim != 2 or A.shape[0] != A.shape[1]:
    raise ValueError("A debe ser una matriz cuadrada.")

if rhs.shape != (A.shape[0],):
    raise ValueError("rhs debe contener un valor por cada ecuación.")

try:
    T = np.linalg.solve(A, rhs)
except np.linalg.LinAlgError as exc:
    raise RuntimeError("El sistema no pudo resolverse de forma única.") from exc

residuos = A @ T - rhs
residuo_max = np.linalg.norm(residuos, ord=np.inf)

for i, valor in enumerate(T, start=1):
    print(f"T{i} = {valor:.4f} °C")

print("Residuos:", residuos)
print(f"Residuo máximo = {residuo_max:.3e}")

Para un sistema lineal se resuelve directamente AT = b. No es necesario calcular primero la inversa de A: además de añadir operaciones, esa ruta puede amplificar errores numéricos.

4.2.2 GNU Octave

GNU Octave permite resolver el mismo sistema con el operador de división matricial izquierda. El residuo vuelve a calcularse después de obtener las temperaturas:

A = [ 4.0  -1.0  -1.0   0.0;
     -2.0   4.5   0.0  -1.0;
     -2.0   0.0   4.2  -1.0;
      0.0  -2.0  -2.0   4.7 ];

rhs = [300.0; 212.5; 106.0; 18.5];

if rows(A) != columns(A)
  error("A debe ser una matriz cuadrada.");
endif

if rows(rhs) != rows(A)
  error("rhs debe contener un valor por cada ecuación.");
endif

T = A \ rhs;
residuos = A*T - rhs;
residuo_max = norm(residuos, Inf);

disp(T)
disp(residuos)
printf("Residuo máximo = %.3e\n", residuo_max)

4.2.3 Qué herramienta conviene utilizar

4.2.4 Cuando la malla crece

En una malla pequeña, almacenar todos los coeficientes en una matriz densa es suficiente. En una malla realista, cada nodo se conecta únicamente con unos pocos vecinos y la mayoría de las posiciones de la matriz son cero. En ese caso conviene utilizar matrices dispersas, conservar una numeración nodal sistemática y automatizar la generación de cada fila a partir del tipo de nodo y de sus condiciones de frontera.

5. Extensiones y consideraciones prácticas

5.1 Extensión al régimen transitorio

En estado estacionario, la energía almacenada en el volumen de control no cambia. En régimen transitorio, la variación de temperatura entre dos instantes debe incorporarse al balance:

El superíndice i representa el instante conocido y i + 1 el nuevo instante que se desea calcular.

Para una malla cuadrada y propiedades constantes se definen la difusividad térmica y el número de Fourier de malla:

El número de Fourier relaciona el avance temporal con la rapidez de difusión térmica a través de la malla.

5.1.1 Método explícito

En la formulación explícita, todos los flujos del lado derecho se evalúan con temperaturas conocidas del instante i. Para un nodo interior de una malla cuadrada, sin generación:

La temperatura nueva se obtiene directamente a partir de valores conocidos del instante anterior.

Su principal ventaja es la sencillez: no exige resolver un sistema simultáneo en cada intervalo. Sin embargo, el paso temporal está limitado por la estabilidad. Para el nodo interior bidimensional anterior debe cumplirse:

Criterio específico para un nodo interior bidimensional, una malla cuadrada y propiedades constantes.

5.1.2 Método implícito

En la formulación implícita, los flujos conductivos se evalúan con las temperaturas desconocidas del instante i + 1. Para el mismo nodo interior:

Las temperaturas del nuevo instante deben resolverse de manera simultánea.

El método implícito no presenta la misma restricción de estabilidad del método explícito, pero eso no significa que cualquier paso temporal produzca una solución precisa. Un Δt demasiado grande puede suavizar o desplazar la evolución térmica y ocultar cambios rápidos.

5.2 Errores frecuentes y cómo evitarlos

La mayoría de los fallos no aparece al resolver el sistema, sino al representar incorrectamente la geometría, omitir un mecanismo o interpretar como válida una salida numérica que no ha sido verificada. La siguiente lista reúne los puntos que conviene revisar antes de aceptar un resultado.

6. Conclusión

El análisis nodal deja de parecer una colección de ecuaciones especiales cuando se conserva una sola lógica: delimitar la región representada por el nodo y aplicar sobre ella conservación de la energía. La ubicación determina si el volumen es completo, medio o un cuarto; las fronteras indican qué mecanismos deben incorporarse; y las propiedades del material establecen el peso de cada conexión térmica.

Temperatura prescrita, flujo impuesto, aislamiento, convección, radiación, generación interna, interfaces y esquinas no son problemas independientes. Son variaciones del mismo balance general. El paso al régimen transitorio tampoco rompe esa estructura: únicamente añade el almacenamiento de energía y obliga a controlar la discretización temporal.

La utilidad del método no termina al construir y resolver el sistema. Una solución numérica adquiere valor ingenieril cuando supera comprobaciones físicas, algebraicas, energéticas y de independencia de malla. Solo entonces puede utilizarse para interpretar el proceso, comparar alternativas o ampliar el modelo a geometrías y mecanismos más complejos.

El tema puede resultar difícil, sobre todo cuando coinciden varias condiciones de frontera. Dominarlo exige práctica, paciencia y volver sobre los balances hasta que cada término tenga sentido físico. Esta guía reúne los fundamentos esenciales, pero no sustituye un tratamiento especializado ni desarrolla con amplitud todas las formulaciones posibles.