Una PDE, tres caminos#

Al poner una PDE en una computadora, lo primero es dónde se sienta la incógnita. La ecuación puede ser la misma. Valor puntual, coeficiente de base o promedio de celda: tres técnicas distintas.

Basta una ecuación de advección-difusión 1D.

$$\frac{\partial u}{\partial t} + a\frac{\partial u}{\partial x} = \nu\frac{\partial^2 u}{\partial x^2}$$

El menú cambia entre FDM, FEM y FVM, y los grados de libertad se mueven entre nodos, coeficientes de base y centros de celda.

FDM alinea los puntos, mientras FVM deja un valor en el centro de cada celda. FEM apoya el coeficiente en el nodo y deja que la base cruce el elemento.

---

1. Método de Diferencias Finitas (Finite Difference Method, FDM)#

Idea central#

“

Aproximar las derivadas directamente como diferencias.

Es el enfoque más intuitivo. Se aproximan las derivadas realizando una expansión de Taylor de los valores de la función en los puntos de la malla (nodos).

Origen matemático: Expansión de Taylor#

Expansión de Taylor en el punto $x_i$:

$$u(x_i + \Delta x) = u(x_i) + \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2}\frac{\partial^2 u}{\partial x^2}\bigg|_i + \cdots$$

A partir de esto, se derivan las aproximaciones por diferencias:

Diferencia hacia adelante (Forward):

$$\frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_{i+1} - u_i}{\Delta x} + O(\Delta x)$$

Diferencia central (Central):

$$\frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_{i+1} - u_{i-1}}{2\Delta x} + O(\Delta x^2)$$

Segunda derivada:

$$\frac{\partial^2 u}{\partial x^2}\bigg|_i \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2} + O(\Delta x^2)$$

Aplicación a la ecuación de advección-difusión#

Usando diferencias centrales:

$$\frac{du_i}{dt} + a\frac{u_{i+1} - u_{i-1}}{2\Delta x} = \nu\frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2}$$

Esto resulta en un sistema de ecuaciones diferenciales ordinarias (ODE) acopladas para las incógnitas $u_i$ en los puntos de la malla.

Ventajas y desventajas#

Ventajas:

• Concepto simple e implementación sencilla.
• Altamente eficiente en mallas estructuradas (structured grids).
• Fácil de lograr alta precisión (esquemas compactos, esquemas tipo espectral).
• Estructura de matriz limpia en mallas ortogonales (matriz banda).

Desventajas:

• Difícil de aplicar a mallas no estructuradas (unstructured grids): esta es una limitación crítica.
• La generación de mallas para geometrías complejas es complicada.
• No satisface automáticamente las leyes de conservación.

Áreas de aplicación representativas#

• DNS/LES (Simulaciones de turbulencia basadas en mallas ortogonales).
• Modelos meteorológicos y oceánicos (mallas estructuradas).
• Simulación de propagación de ondas sísmicas.
• Cálculos de alta precisión basados en esquemas compactos.

---

2. Método de Elementos Finitos (Finite Element Method, FEM)#

Idea central#

“

Aproximar la solución como una combinación lineal de funciones de base y minimizar el residuo pesado.

El FEM no resuelve la ecuación diferencial directamente. En su lugar, la transforma en una forma integral denominada formulación débil (weak formulation).

Origen matemático: Forma débil (Weak Form)#

Se multiplica la PDE original (forma fuerte o strong form) por una función de prueba (test function) $w$ y se integra:

$$\int_\Omega w \left(\frac{\partial u}{\partial t} + a\frac{\partial u}{\partial x} - \nu\frac{\partial^2 u}{\partial x^2}\right) dx = 0$$

Aplicando integración por partes al término de difusión:

$$\int_\Omega w \frac{\partial u}{\partial t}\,dx + \int_\Omega w\, a\frac{\partial u}{\partial x}\,dx + \int_\Omega \nu\frac{\partial w}{\partial x}\frac{\partial u}{\partial x}\,dx = \left[\nu w \frac{\partial u}{\partial x}\right]_\Gamma$$

Esta es la forma débil. Los cambios que importan son estos:

• Se requería la segunda derivada original, pero tras la integración por partes, solo se necesita la primera derivada.
• Se relajan los requisitos de continuidad de la solución ($C^1 \to C^0$).
• Las condiciones de contorno se incluyen de forma natural (BC de Neumann = lado derecho).

Aproximación de Galerkin#

Se expresa la solución como una combinación lineal de funciones de base $\phi_j$:

$$u^h(x, t) = \sum_{j=1}^{N} U_j(t)\,\phi_j(x)$$

En el método de Galerkin, las funciones de prueba y las funciones de base se eligen del mismo espacio ($w = \phi_i$):

$$\sum_j \left(\int_\Omega \phi_i \phi_j\,dx\right) \frac{dU_j}{dt} + \sum_j \left(\int_\Omega \phi_i\, a\, \phi_j'\,dx + \int_\Omega \nu\, \phi_i'\, \phi_j'\,dx\right) U_j = 0$$

En forma matricial:

$$\mathbf{M}\frac{d\mathbf{U}}{dt} + \mathbf{K}\mathbf{U} = \mathbf{f}$$

• $\mathbf{M}$: Matriz de masa (mass matrix).
• $\mathbf{K}$: Matriz de rigidez (stiffness matrix).
• $\mathbf{f}$: Vector de contorno/fuente.

Elección de funciones de base#

Lo más habitual es utilizar elementos basados en polinomios de Lagrange:

Ventajas y desventajas#

Ventajas:

• Natural para mallas no estructuradas: maneja formas complejas con mallas triangulares o tetraédricas.
• Es posible realizar estimaciones de error matemáticamente rigurosas (a priori / a posteriori).
• Se combina bien con el refinamiento de malla adaptativo.
• Flexibilidad para elegir entre refinamiento p (aumentar el grado del polinomio) y refinamiento h (refinamiento de malla).

Desventajas:

• No satisface las leyes de conservación directamente (Galerkin estándar).
• Inestable en problemas dominados por la advección (oscilación): requiere estabilización como SUPG o GLS.
• Requiere la inversión de la matriz de masa (o mass lumping).
• Mayor complejidad de implementación en comparación con FDM/FVM.

Áreas de aplicación representativas#

• Mecánica de sólidos / Análisis estructural (hogar original del FEM).
• Análisis de campos electromagnéticos.
• Transferencia de calor.
• Biomecánica.
• Geomecánica.

---

3. Método de Volúmenes Finitos (Finite Volume Method, FVM)#

Idea central#

“

Integrar las leyes de conservación sobre volúmenes de control y equilibrar los flujos que pasan a través de los límites del volumen.

El FVM refleja las leyes de conservación física de la forma más directa.

Origen matemático: Ley de conservación integral#

Se escribe la ecuación de advección-difusión en forma integral:

$$\frac{d}{dt}\int_{V_i} u\,dx + \oint_{\partial V_i} \left(au - \nu\frac{\partial u}{\partial x}\right) \cdot n\, dS = 0$$

Para un volumen de control $V_i = [x_{i-1/2},\, x_{i+1/2}]$:

$$\frac{d\bar{u}_i}{dt} = -\frac{1}{\Delta x}\left(\hat{F}_{i+1/2} - \hat{F}_{i-1/2}\right)$$

Donde $\bar{u}_i$ es el valor promedio de la celda y $\hat{F}_{i+1/2}$ es el flujo numérico en el límite de la celda.

Determinación de los flujos numéricos#

El núcleo del FVM es cómo calcular los flujos en los límites de las celdas:

Término de advección: basado en simuladores de Riemann o esquemas upwind:

$$\hat{F}^{adv}_{i+1/2} = \text{RiemannSolver}(\bar{u}_i,\, \bar{u}_{i+1})$$

Término de difusión: diferencia central:

$$\hat{F}^{diff}_{i+1/2} = -\nu\frac{\bar{u}_{i+1} - \bar{u}_i}{\Delta x}$$

Reconstrucción#

Los promedios de celda por sí solos ofrecen solo precisión de primer orden. Para una mayor precisión, se debe reconstruir la distribución dentro de la celda:

MUSCL (2.º orden):

$$u_{i+1/2}^L = \bar{u}_i + \frac{1}{2}\psi(r_i)\,(\bar{u}_i - \bar{u}_{i-1})$$

Donde $\psi$ es un limitador de pendiente y $r_i$ es la relación de gradientes sucesivos.

WENO (5.º orden):

Logra una precisión de quinto orden mediante el promedio ponderado no lineal de tres candidatos a polinomios de segundo orden:

$$u_{i+1/2}^L = \sum_{k=0}^{2} \omega_k\, q_k(x_{i+1/2})$$

Los pesos $\omega_k$ se determinan según los indicadores de suavidad $\beta_k$, otorgando automáticamente mayor peso a los esquemas suaves cerca de las discontinuidades.

La difusión de cara escrita arriba es $-\nu(\bar{u}_{i+1}-\bar{u}_i)/\Delta x$, un espaciado ortogonal. Esa diferencia es la derivada normal solo cuando la normal sigue la línea de la malla. Una cara torcida no saca la dirección de los promedios vecinos nada más. La corrección de difusión no ortogonal cierra ese hueco. Sin ella, una actualización conservativa puede empujar la difusión fuera de eje.

La diferencia hacia adelante llevaba $O(\Delta x)$. La central y el estencil de la segunda derivada llevaban $O(\Delta x^2)$. La hacia adelante conserva el término de Taylor $\Delta x^2/2$. La central resta $u_{i+1}$ y $u_{i-1}$ y lo cancela. Al partir el espaciado a la mitad, el error hacia adelante baja cerca de la mitad. Las piezas $O(\Delta x^2)$ bajan cerca de un cuarto. Promedios de celda solos, sin reconstrucción, se quedan en primer orden. El orden extra se compra al subir hasta la cara.

El $1/2$ de MUSCL es el factor que levanta el promedio de celda, en línea recta, hasta la cara. Con $\psi(r_i)=1$ se guarda toda la pendiente. Con $\psi=0$ el valor de cara vuelve a $\bar{u}_i$ y el flujo es otra vez de primer orden. WENO suma solo tres candidatos, $k=0,1,2$. Cada $q_k$ es una cuadrática. $\omega_k$ recorta peso donde $\beta_k$ es grande. "Quinto orden en un tramo suave" es el orden de esa suma pesada. Si $\psi$ cae a 0, la fórmula de quinto orden sigue escrita, pero esa cara bajó a primer orden. Un triángulo P1 tiene 3 nodos, P2 tiene 6 y P3 tiene 10. Subir el grado de 1 a 3 empuja las incógnitas por elemento de 3 hasta pasar de 10. En la forma débil, $\mathbf{M}$ guarda $\int \phi_i \phi_j\,dx$ y $\mathbf{K}$ junta $\phi_i \phi_j'$ con $\nu\,\phi_i'\phi_j'$. Con $\nu$ positivo el bloque de difusión queda casi simétrico. La velocidad advectiva $a$ deja la parte antisimétrica.

Ventajas y desventajas#

Ventajas:

• Las leyes de conservación se satisfacen exactamente a nivel discreto: una ventaja decisiva en CFD.
• Excelente para capturar ondas de choque y discontinuidades.
• Aplicable a mallas no estructuradas.
• El pensamiento basado en flujos se alinea directamente con la intuición física.

Desventajas:

• Lograr una alta precisión es más difícil que en el FEM (la reconstrucción se vuelve compleja).
• Puede ser menos eficiente que el FEM para ecuaciones de difusión puras.
• La implementación de la reconstrucción de alto orden en mallas no estructuradas es complicada.

Áreas de aplicación representativas#

• Toda el área de CFD (OpenFOAM, ANSYS Fluent, SU2).
• Flujo compresible/incompresible.
• Flujo multifásico, combustión, flujos reactivos.
• Parte de los núcleos dinámicos en modelos meteorológicos y climáticos.

---

Comparación clave#

Puntos de partida matemáticos#

¿"Qué" se discretiza?#

Esta es la diferencia más fundamental:

• FDM: Discretiza el operador diferencial ($\partial/\partial x$).
• FEM: Discretiza el espacio de la solución.
• FVM: Discretiza el flujo en las ecuaciones integrales.

Ubicación de las incógnitas#

• FDM: Valores puntuales en los puntos de la malla (nodos).
• FEM: Valores puntuales en los nodos (coeficientes de las funciones de base).
• FVM: Valores promedio de celda.

Comparación integral#

---

Técnicas modernas: difuminando las fronteras#

Recientemente, se han investigado activamente métodos híbridos que combinan las ventajas de las tres técnicas:

Galerkin discontinuo (Discontinuous Galerkin, DG)#

Combina las funciones de base del FEM con el concepto de flujo del FVM. Aproxima la solución con polinomios dentro de cada celda (FEM) y calcula los flujos en los límites de las celdas mediante simuladores de Riemann (FVM).

$$\int_{V_i} w \frac{\partial u^h}{\partial t}\,dx - \int_{V_i} \frac{\partial w}{\partial x} F(u^h)\,dx + \hat{F}_{i+1/2} w_{i+1/2}^- - \hat{F}_{i-1/2} w_{i-1/2}^+ = 0$$

Satisface simultáneamente la conservación, la precisión de alto orden y la capacidad de mallas no estructuradas. No obstante, los costos computacionales son elevados y se requieren limitadores cerca de los choques.

Diferencia espectral / Reconstrucción de flujo#

Similar al DG pero se resuelve en forma diferencial sin integración para aumentar la eficiencia. Se combina muy bien con la aceleración por GPU y está ganando atención para los simuladores de CFD de próxima generación.

Métodos sin malla (Meshless Methods como SPH)#

Un enfoque que elimina la malla por completo, aproximando basándose en partículas. Fuerte para flujos de superficie libre y problemas de grandes deformaciones, pero existen problemas de precisión y consistencia.

---

Si la integral tiene que volver entera, se deja el estencil puro#

Un DNS en una caja, con malla ortogonal, vive con el estencil central $O(\Delta x^2)$. Tres nodos arman la segunda derivada y la matriz sigue en banda. Si la geometría se dobla y la malla deja de ser estructurada, ese estencil ya no nombra un vecino. La forma débil se traga la geometría. Aun así, $\mathbf{M}$ y $\mathbf{K}$ del Galerkin estándar no prometen que los flujos se cancelen celda a celda.

Cuando pasa un choque y la integral tiene que regresar como el mismo número, el esquema es FVM. El flujo de cara $\hat{F}_{i+1/2}$ es pérdida para la celda izquierda y ganancia para la derecha, así que las caras interiores se anulan en la suma global. Orden alto y conservación en un solo esquema significan un polinomio dentro de la celda y un flujo de Riemann en la cara, como en una discretización de Galerkin discontinuo. Sin limitador, ese híbrido también oscila delante del choque.