Lectura crítica de simulaciones

Criterios de validación física, análisis de sensibilidad, conservación de magnitudes y detección de artefactos.

1. El Paradigma de Verificación y Validación (V&V)

En la ciencia e ingeniería computacional contemporánea, una simulación no constituye una prueba empírica por sí misma, sino la resolución numérica de un modelo matemático idealizado. Para establecer la credibilidad y exactitud científica de una simulación física, se recurre a la disciplina formal de Verificación y Validación (V&V):

A. Verificación Numérica: "¿Se están resolviendo las ecuaciones correctamente?"

La verificación es un proceso estrictamente matemático y computacional. Determina si los algoritmos programados resuelven el modelo matemático con la precisión esperada, libre de errores de código (bugs), errores de discretización espaciotemporal y desbordamientos de redondeo en punto flotante. Comprende dos etapas:

1. Verificación del código: Comprobación del orden asintótico teórico de convergencia mediante el Método de Soluciones Manufacturadas (MMS) o soluciones exactas.
2. Verificación de la solución: Estimación cuantitativa de la incertidumbre numérica de la solución discreta respecto a la solución analítica continua del mismo modelo.

B. Validación Física: "¿Se están resolviendo las ecuaciones correctas?"

La validación es un proceso físico y empírico. Determina en qué medida el modelo matemático subyacente y sus hipótesis constitutivas (por ejemplo, fluido incompresible, gas ideal, linealidad elástica, potencial de Coulomb puntual) reproducen fielmente los fenómenos observados en experimentos reales dentro del rango de operación previsto.

$\text{Mundo Real (Experimento)} \xrightarrow{\text{Validación}} \text{Modelo Conceptual/Matemático} \xrightarrow{\text{Verificación}} \text{Simulación Numérica}$

2. Metodología de Convergencia de Malla e Índice GCI de Roache

Para demostrar que una solución discreta es independiente de la resolución espacial o temporal, se debe realizar un estudio de convergencia sistemático empleando un conjunto de al menos tres mallas refinadas sucesivamente con un factor de escala constante $r = h_{\text{gruesa}} / h_{\text{fina}} > 1.3$.

Sean $f_1, f_2, f_3$ los valores calculados de una magnitud de interés (por ejemplo, energía total, fuerza de arrastre o presión en un punto crítico) en la malla fina ($h_1$), intermedia ($h_2 = r h_1$) y gruesa ($h_3 = r^2 h_1$).

Estimación del Orden de Convergencia Aparente ($p$):

$p = \frac{\ln\left(\frac{f_3 - f_2}{f_2 - f_1}\right)}{\ln(r)}$

• Si $p$ coincide con el orden formal teórico del algoritmo discretizado (por ejemplo, $p \approx 2$ para Crank-Nicolson o diferencias centradas), la simulación se halla en el régimen de convergencia asintótica.
• Si $p < 1$ o varía erráticamente, la malla aún es demasiado gruesa para capturar las singularidades o capas límite de la física del problema.

Extrapolación de Richardson:

La extrapolación predice la solución en el límite continuo ideal ($h \to 0$): $f_{\text{extrap}} = f_1 + \frac{f_1 - f_2}{r^p - 1}$

Índice de Convergencia de Malla (Grid Convergence Index - GCI):

Estandarizado por la ASME y AIAA, el GCI provee una banda de incertidumbre porcentual con un factor de seguridad $F_s = 1.25$ para tres mallas:

$\text{GCI}_{12} = \frac{F_s \, |e_a|}{r^p - 1} \times 100\%, \quad e_a = \frac{f_1 - f_2}{f_1}$

3. Leyes de Conservación e Invariantes Físicos como Diagnóstico

Un código numérico puede converger formalmente pero violar principios físicos inviolables debido a simplificaciones algebraicas o integraciones asimétricas:

1. Conservación de la Energía: En sistemas conservativos no forzados, la deriva relativa $\Delta E(t) / E_0$ debe estar estrictamente acotada o converger a cero cuando $h \to 0$.
2. Conservación de Cantidad de Movimiento e Invariancia Galileana: Fuerzas internas no pueden acelerar el centro de masas del sistema.
3. Unitariedad Cuántica y Conservación de Probabilidad: En mecánica cuántica, $\int |\psi(\vec{x}, t)|^2 d^3x \equiv 1$.
4. Conservación de Masa y Continuidad: El flujo neto en los contornos debe igualar la tasa de acumulación interna $\oint \vec{J} \cdot d\vec{A} = -\frac{d}{dt} \int \rho \, dV$.

4. Detección y Caracterización de Artefactos Numéricos

Al examinar resultados de simulaciones complejas, el físico debe identificar y aislar artefactos espurios originados exclusivamente por la computación:

5. Cuantificación de Incertidumbre (UQ): Aleatoria vs Epistémica

La lectura crítica distingue dos fuentes de incertidumbre que limitan las conclusiones de una simulación:

• Incertidumbre Aleatoria (Irreducible): Variabilidad estocástica inherente al sistema físico (ruido térmico Browniano, turbulencia plenamente desarrollada, fluctuaciones cuánticas). Se cuantifica mediante muestreo de Monte Carlo y funciones de densidad de probabilidad.
• Incertidumbre Epistémica (Reducible): Falta de información o ignorancia sobre los parámetros del modelo (conductividad térmica exacta del material, rugosidad superficial no medida, coeficientes de reacción química estimados). Se reduce con experimentos específicos adicionales.

Protocolo de Auditoría y Lectura Crítica de Simulaciones Físicas

Metodología Sistemática de 5 Fases

Para evaluar con rigor científico los resultados de una simulación computacional (ya sea propia, de la literatura o de software comercial), se debe ejecutar la siguiente auditoría metódica:

Fase 1: Auditoría del Modelo Conceptual y Rango de Validez

1. Identificar las simplificaciones teóricas del modelo: ¿Es el flujo laminar o turbulento? ¿Se asume linealidad elástica en presencia de deformaciones finitas? ¿Se desprecian efectos relativistas o cuánticos válidamente?
2. Verificar si los números adimensionales del sistema (Reynolds $\text{Re}$, Mach $\text{Ma}$, Knudsen $\text{Kn}$, Péclet $\text{Pe}$) se encuentran dentro del régimen en el que las ecuaciones gobernantes son físicamente válidas.

Fase 2: Comprobación de Convergencia de Malla e Independencia Espaciotemporal

1. No aceptar jamás una simulación realizada en una única malla o con un único paso de tiempo.
2. Exigir un estudio de al menos tres mallas ($h_1 < h_2 < h_3$ con $r \ge 1.3$).
3. Calcular el orden aparente $p$ y el índice GCI de Roache:
   • Si $\text{GCI}_{12} < 1\% - 2\%$, la incertidumbre por discretización está adecuadamente controlada.
   • Si $\text{GCI}_{12} > 10\%$, los resultados cuantitativos no tienen validez conclusiva.

Fase 3: Telemetría de Leyes de Conservación e Invariantes

1. Graficar la evolución temporal de los invariantes fundamentales: $\Delta E(t) = \frac{|E(t) - E(0)|}{|E(0)|}, \quad \Delta M(t) = \frac{|M(t) - M(0)|}{|M(0)|}$
2. Verificar el balance global de flujos en las fronteras: el calor/masa que entra menos el que sale debe coincidir numéricamente con la tasa de variación temporal del volumen interior.

Fase 4: Búsqueda Activa de Artefactos Numéricos

1. Inspección de simetrías: Si la geometría y las condiciones de contorno son simétricas, ¿presenta la solución asimetrías sospechosas inducidas por mallas no estructuradas o ruido numérico?
2. Inspección de bordes artificiales: Verificar si ondas salientes o vórtices rebotan espuriamente en los límites numéricos de contorno (falta de condiciones absorbentes tipo PML o no reflectivas).
3. Comprobación de suavidad: Buscar oscilaciones espurias de alta frecuencia (efecto tablero de ajedrez en presión-velocidad, o rizado de Gibbs cerca de choques).

Fase 5: Validación Cruzada Frente a Benchmarks y Datos Experimentales

1. Contrastar el código frente a casos canónicos con solución analítica exacta (onda plana acústica, flujo de Poiseuille, oscilador armónico simple).
2. Comparar cuantitativamente los resultados simulados con mediciones experimentales independientes, incluyendo las barras de error empíricas y las bandas de incertidumbre GCI.

Ejemplos Resueltos: Auditoría Cuantitativa de Simulaciones

Ejemplo 1: Cálculo del Índice de Convergencia de Malla (GCI) para el Coeficiente de Arrastre

Un equipo de aerodinámica simula el flujo de aire alrededor de un perfil alar para determinar el coeficiente de sustentación $C_L$. Ejecutan tres simulaciones en mallas sucesivamente refinadas con factor $r = \sqrt{2} \approx 1.414$:

• Malla 3 (gruesa, $N_3 = 25000$ nodos): $C_{L,3} = 1.0820$
• Malla 2 (media, $N_2 = 50000$ nodos): $C_{L,2} = 1.1240$
• Malla 1 (fina, $N_1 = 100000$ nodos): $C_{L,1} = 1.1350$

1. Determinación de los errores relativos aparentes:

$e_{23} = \frac{C_{L,2} - C_{L,3}}{C_{L,2}} = \frac{1.1240 - 1.0820}{1.1240} = \frac{0.0420}{1.1240} = 0.03737 \, (3.74\%)$ $e_{12} = \frac{C_{L,1} - C_{L,2}}{C_{L,1}} = \frac{1.1350 - 1.1240}{1.1350} = \frac{0.0110}{1.1350} = 0.00969 \, (0.97\%)$

2. Cálculo del orden de convergencia aparente ($p$):

$p = \frac{\ln\left(\frac{C_{L,3} - C_{L,2}}{C_{L,2} - C_{L,1}}\right)}{\ln(r)} = \frac{\ln\left(\frac{-0.0420}{-0.0110}\right)}{\ln(1.4142)} = \frac{\ln(3.81818)}{0.34657} = \frac{1.33977}{0.34657} \approx 3.866$

Como $p$ es positivo y los saltos se reducen monótonamente ($f_1 > f_2 > f_3$), la solución está convergiendo de forma monótona hacia un valor límite.

3. Estimación por Extrapolación de Richardson ($h \to 0$):

$C_{L,\text{extrap}} = C_{L,1} + \frac{C_{L,1} - C_{L,2}}{r^p - 1} = 1.1350 + \frac{0.0110}{(1.4142)^{3.866} - 1} = 1.1350 + \frac{0.0110}{3.818 - 1} = 1.1350 + \frac{0.0110}{2.818} = 1.1389$

4. Cálculo del GCI de la malla fina (con factor de seguridad $F_s = 1.25$):

$\text{GCI}_{12} = \frac{F_s \, |e_{12}|}{r^p - 1} = \frac{1.25 \times 0.00969}{2.818} = \frac{0.01211}{2.818} = 0.00430 \implies 0.43\%$

Conclusión Científica: El valor simulado en la malla fina es $C_L = 1.1350 \pm 0.0049$ (incertidumbre numérica de discretización del $0.43\%$). La simulación se encuentra sólidamente verificada.

Ejemplo 2: Diagnóstico de Fricción Física vs Disipación Numérica Espuria

Se simula un oscilador armónico simple ideal de masa $m = 1\text{ kg}$ y resorte $k = 100\text{ N/m}$ (frecuencia natural $\omega = 10\text{ rad/s}$, período $T = 0.628\text{ s}$). No hay amortiguamiento físico alguno ($\gamma = 0$). Energía inicial: $E(0) = 50.0\text{ J}$.

Al cabo de $t = 10\text{ segundos}$ (aproximadamente 16 períodos), el archivo de resultados registra:

• Con integrador A: $E(10) = 42.1\text{ J}$
• Con integrador B: $E(10) = 50.000003\text{ J}$
• Con integrador C: $E(10) = 118.5\text{ J}$

Diagnóstico Físico-Numérico:

1. Integrador A (Euler Implícito o Disipación Numérica Severa): Ha disipado $7.9\text{ J}$ ($15.8\%$ de la energía mecánica total) en ausencia total de amortiguadores. El código introduce una fricción artificial que simula falsamente un medio viscoso viscoso inexistente. Diagnóstico: Rechazar para dinámica de sistemas conservativos.
2. Integrador B (Verlet de Velocidades o RK4 con paso adecuado): La energía se conserva con una fluctuación relativa de $6 \times 10^{-8}$. Diagnóstico: Físicamente fiel y matemáticamente consistente.
3. Integrador C (Euler Explícito): La energía aumentó de $50\text{ J}$ a $118.5\text{ J}$ ($+137\%$). El algoritmo inyecta energía no física debido a que el determinante de su matriz de transferencia excede la unidad. Diagnóstico: Completamente espurio y desestabilizante.

Errores Comunes en la Lectura y Evaluación de Simulaciones

1. Confundir Atractivo Visual con Precisión Científica

• Error: Concluir que una simulación 3D con renderizado fotorrealista y colores vistosos es físicamente correcta.
• Consecuencia: Se aceptan como válidos modelos que pueden contener errores flagrantes de discretización, inestabilidades enmascaradas o violaciones termodinámicas severas.
• Corrección: Exigir siempre gráficos cuantitativos de convergencia, telemetría de energía y balances de invariantes antes de mirar las animaciones visuales.

2. Realizar Ajuste Paramétrico Espurio (Overfitting Numérico)

• Error: Modificar arbitrariamente coeficientes empíricos no medidos (fricciones, viscosidades turbulentas) hasta que la curva simulada encaje con los puntos experimentales.
• Consecuencia: La simulación pierde toda capacidad predictiva; cualquier cambio en las condiciones operativas destruirá la concordancia.
• Corrección: Separar estrictamente la calibración de parámetros de la validación experimental independiente en condiciones no vistas.

3. Aceptar una Simulación Basada en una Sola Malla

• Error: Publicar o dar por buenos resultados numéricos obtenidos en una única resolución de malla sin verificar independencia espacial.
• Consecuencia: Imposibilidad de saber si los vórtices o fluctuaciones observadas son físicas o producto del tamaño de la celda elegida.
• Corrección: Realizar siempre estudios de convergencia con al menos tres niveles de malla y reportar el índice GCI.

4. Confundir Convergencia Iterativa del Solver con Convergencia Física

• Error: Pensar que porque los residuos algebraicos cayeron por debajo de $10^{-6}$, el resultado físico es exacto con 6 cifras decimales.
• Consecuencia: Los residuos evalúan la resolución del sistema algebraico discreto, no el error de aproximación entre la malla discreta y la realidad física continua.
• Corrección: Medir el error de discretización espacial y temporal mediante refinamiento sistemático.

5. Ignorar la Reflexión en Contornos Computacionales Artificiales

• Error: Truncar el dominio de cálculo cerca de la zona de interés sin implementar condiciones de frontera no reflectivas o capas absorbentes (PML).
• Consecuencia: Ondas acústicas, electromagnéticas o elásticas rebotan en los límites de la caja de simulación y regresan a la zona de estudio, falseando resonancias.
• Corrección: Alejar las fronteras artificiales y verificar la absorción perfecta de ondas salientes.

Aplicaciones Físicas de la Auditoría y Lectura Crítica de Simulaciones

1. Certificación de Seguridad Aerodinámica y Aeroespacial

En la industria aeroespacial (NASA, ESA, Boeing, Airbus), la certificación de aeronaves y cohetes exige protocolos estrictos de V&V: $\text{GCI} < 1.0\% \text{ en coeficientes de sustentación } C_L \text{ y arrastre } C_D$ Antes de realizar vuelos de prueba con tripulación, las simulaciones de entrada atmosférica y choque hipersónico deben demostrar que la temperatura de la pantalla térmica no presenta artefactos de discretización que puedan subestimar la ablación del material.

2. Modelización Climática Global y Predicción Meteorológica

En los Modelos de Circulación General (GCMs) que proyectan el cambio climático global:

• Se evalúa críticamente la sensibilidad del clima ante duplicaciones de $\text{CO}_2$.
• La lectura crítica permite separar la variabilidad climática natural (incertidumbre aleatoria) de las simplificaciones en las parametrizaciones de nubes y albedo del hielo (incertidumbre epistémica).

3. Astrofísica de Ondas Gravitacionales y Fusión de Agujeros Negros

En el análisis de señales detectadas por LIGO y Virgo, los patrones de ondas gravitacionales emitidos durante la coalescencia de agujeros negros binarios se generan mediante relatividad numérica: $G_{\mu\nu} = \frac{8\pi G}{c^4} T_{\mu\nu}$ La lectura crítica audita que la energía y el momento angular transportados por la radiación gravitacional cuadrupolar coincidan exactamente con la disminución de masa ADM del sistema binario.

4. Termohidráulica y Seguridad de Reactores Nucleares

En el diseño y licenciamiento de reactores nucleares de fisión y fusión (ITER):

• Se exige la validación cruzada de códigos termohidráulicos (RELAP5, TRACE) frente a incidentes de pérdida de refrigerante (LOCA).
• La auditoría crítica garantiza que la predicción de temperatura pico en las vainas de combustible no esté subestimada por difusión numérica espuria.

5. Biomecánica Cardiovascular y Flujo en Aneurismas

En la planificación quirúrgica de aneurismas cerebrales y colocación de stents vasculares:

• Las simulaciones hemodinámicas calculan la tensión de cizalla en la pared endotelial (WSS).
• La verificación de malla es crítica: un mallado insuficiente en la pared arterial puede subestimar las fuerzas cortantes hasta en un $50\%$, llevando a diagnósticos clínicos erróneos sobre el riesgo de ruptura del vaso.

Problemas Prácticos: Verificación, Validación y Auditoría de Simulaciones

Problema 1: Detección de Falso Orden Asintótico en una Malla No Resuelta

Un estudiante modela la propagación de calor con un esquema que teóricamente es de segundo orden espacial ($\mathcal{O}(\Delta x^2)$). Realiza tres simulaciones con refinamiento $r = 2$:

• Malla gruesa ($h_3 = 0.4\text{ m}$): $T_3 = 320.0\text{ K}$
• Malla media ($h_2 = 0.2\text{ m}$): $T_2 = 312.0\text{ K}$
• Malla fina ($h_1 = 0.1\text{ m}$): $T_1 = 309.0\text{ K}$

1. Calcule el orden aparente de convergencia $p$.
2. Compare $p$ con el orden teórico esperado ($p = 2$). ¿Qué revela esta discrepancia sobre la calidad de la discretización?

Solución Paso a Paso:

1. Calculamos las diferencias entre soluciones consecutivas: $f_3 - f_2 = 320.0 - 312.0 = 8.0\text{ K}$ $f_2 - f_1 = 312.0 - 309.0 = 3.0\text{ K}$ Con ratio de refinamiento $r = 2$: $p = \frac{\ln\left(\frac{8.0}{3.0}\right)}{\ln(2)} = \frac{\ln(2.6667)}{0.69315} = \frac{0.98083}{0.69315} \approx 1.415$
2. El orden formal teórico es $p_{\text{teórico}} = 2.0$, pero el orden observado es apenas $p \approx 1.41$. Esto evidencia que la malla aún no se encuentra en el rango asintótico. Probablemente existen gradientes locales intensos o capas límite térmicas que requieren celdas más finas que $h = 0.1\text{ m}$ para que las derivadas de orden superior de Taylor sean verdaderamente despreciables.

Problema 2: Auditoría de Conservación de Masa en CFD

En una simulación bidimensional de un canal convergente para un líquido incompresible ($\rho = 1000\text{ kg/m}^3$), el caudal másico medido en la entrada es $\dot{m}_{\text{in}} = 15.42\text{ kg/s}$, mientras que en la salida se registra $\dot{m}_{\text{out}} = 14.88\text{ kg/s}$.

1. Calcule el desbalance relativo de masa porcentual.
2. Si el fabricante del software comercial afirma que "la simulación convergió con residuos inferiores a $10^{-4}$", explique por qué dicho criterio de residuos iterativos no garantiza la conservación física de masa.

Solución Paso a Paso:

1. Desbalance relativo: $\Delta \dot{m} = \frac{|\dot{m}_{\text{in}} - \dot{m}_{\text{out}}|}{\dot{m}_{\text{in}}} = \frac{15.42 - 14.88}{15.42} = \frac{0.54}{15.42} \approx 0.0350 \implies 3.50\%$ Se está "destruyendo" artificialmente el $3.5\%$ del fluido dentro del canal.
2. Los residuos algebraicos de un solver evalúan únicamente la discrepancia interna del sistema lineal iterativo ($||A x - b|| < 10^{-4}$), no la consistencia global integral. Si la formulación numérica no es estrictamente conservativa (método de volúmenes finitos con balances exactos en las caras de las celdas), o si la tolerancia de presión-velocidad no es lo bastante estricta, pueden existir pérdidas continuas de masa a pesar de que el solver declare convergencia algebraica exitosa.

Problema 3: Detección de Dispersión Numérica en Ondas Acústicas

En una simulación acústica donde viaja un paquete de ondas de pulso gaussiano centrado en $f = 1000\text{ Hz}$, se observa que tras recorrer $20\text{ m}$ aparecen oscilaciones sinusoidales detrás del paquete que no estaban presentes al inicio.

1. Identifique el fenómeno numérico causante de estas oscilaciones parásitas.
2. Indique dos medidas correctivas para eliminar este artefacto sin alterar la física del medio.

Solución Paso a Paso:

1. El fenómeno es la dispersión numérica. Ocurre porque la velocidad de fase numérica en una malla discreta depende de la longitud de onda: las componentes de alta frecuencia viajan más despacio que las de baja frecuencia cuando el número de Courant es $C < 1$ o cuando hay menos de 10-15 nodos por longitud de onda ($\lambda / \Delta x < 15$). Las componentes rezagadas forman una estela oscilatoria parásita (ringing).
2. Medidas correctivas:
   • Refinar la malla espacial para garantizar al menos 20 nodos por longitud de onda mínima ($\Delta x \le \lambda_{\min} / 20$).
   • Ajustar el paso temporal para que el número de Courant esté más próximo a 1 (donde la dispersión se minimiza en mallas uniformes), o emplear esquemas de orden superior (como diferencias finitas compactas o espectrales).