(CFD) a menudo se ve como una caja negra de software comercial complejo. Sin embargo, implementar un solucionador “desde cero” es una de las formas más poderosas de aprender la física del movimiento de fluidos. Comencé esto como un proyecto personal y, como parte de un curso de Biofísica, lo aproveché como una oportunidad para comprender finalmente cómo funcionan estas hermosas simulaciones.
Esta guía está diseñada para científicos e ingenieros de datos que desean ir más allá de las bibliotecas de alto nivel y comprender la mecánica subyacente de las simulaciones numéricas mediante la traducción de ecuaciones diferenciales parciales a código Python discretizado. También exploraremos conceptos fundamentales de programación como operaciones vectorizadas con NumPy y convergencia estocástica, que son habilidades esenciales para todos los interesados en arquitecturas de aprendizaje automático y computación científica más amplias.
Analizaremos la derivación y la implementación en Python de un solucionador simple e incompresible de Navier-Stokes (NS). Y luego aplicaremos este solucionador para simular el flujo de aire alrededor del perfil del ala de un pájaro.
La física: Navier-Stokes incompresible
Las ecuaciones fundamentales de CFD son las ecuaciones de Navier-Stokes, que describen cómo evolucionan la velocidad y la presión en un fluido. Para un vuelo estable (como un pájaro planeando), asumimos que el aire es incompresible (densidad constante) y laminar. Pueden entenderse simplemente como la ley del movimiento de Newton, pero para un elemento infinitesimal de fluido, con las fuerzas que lo afectan. Esto es principalmente presión y viscosidad, pero dependiendo del contexto, podría agregar gravedad, tensiones mecánicas o incluso electromagnetismo si lo desea. El autor puede dar fe de que esto no es recomendable para un primer proyecto.
Las ecuaciones en forma vectorial son:
\[
\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{v} \\
\nabla \cdot \mathbf{v} = 0
\]
Dónde:
v: Campo de velocidades (u,v) p: Presión ρ: Densidad del fluido ν: Viscosidad cinemática
La primera ecuación (Momentum) equilibra la inercia con los gradientes de presión y la difusión viscosa. La segunda ecuación (Continuidad) obliga a que la densidad del fluido permanezca constante.
El problema del acoplamiento de presión
Un desafío importante en CFD es que la presión y la velocidad están acopladas: el campo de presión debe ajustarse constantemente para garantizar que el fluido permanezca incompresible.
Para resolver esto, derivamos una ecuación de presión-Poisson tomando la divergencia de la ecuación del momento. En un solucionador discretizado, resolvemos esta ecuación de Poisson en cada paso de tiempo para actualizar la presión, asegurando que el campo de velocidad permanezca libre de divergencias.
Discretización: de las matemáticas a la cuadrícula
Para resolver estas ecuaciones en una computadora, utilizamos esquemas de diferencias finitas en una cuadrícula uniforme.
Hora: Diferencia hacia delante (Euler explícito). Advección (términos no lineales): diferencia hacia atrás/contra el viento (para estabilidad). Difusión y presión: diferencia central.
Por ejemplo, la fórmula de actualización para el componente u (velocidad x) se ve así en forma de diferencias finitas
\[
u_{i,j}^n \frac{u_{i,j}^n - u_{i-1,j}^n}{\Delta x}
\]
En código, el término de advección u∂x∂u usa una diferencia hacia atrás:
\[
u_{i,j}^n \frac{u_{i,j}^n - u_{i-1,j}^n}{\Delta x}
\]
La implementación de Python
La implementación se realiza en cuatro pasos distintos utilizando matrices NumPy.
1. Inicialización
Definimos el tamaño de la cuadrícula (nx, ny), el paso de tiempo (dt) y los parámetros físicos (rho, nu). Inicializamos los campos de velocidad (u,v) y presión (p) a ceros o un flujo uniforme.
2. La geometría del ala (límite inmerso)
Para simular un ala en una cuadrícula cartesiana, necesitamos marcar qué puntos de la cuadrícula se encuentran dentro del ala sólida.
Cargamos una malla de ala (por ejemplo, desde un archivo STL). Creamos una matriz de máscara booleana donde Verdadero indica un punto dentro del ala. Durante la simulación, forzamos la velocidad a cero en estos puntos enmascarados (condición de no deslizamiento/no penetración).
3. El bucle principal del solucionador
El bucle central se repite hasta que la solución alcanza un estado estable. Los pasos son:
Construya el término fuente (b): calcule la divergencia de los términos de velocidad. Resolver presión: Resuelva la ecuación de Poisson para p usando la iteración de Jacobi. Actualizar velocidad: utilice la nueva presión para actualizar u y v. Aplicar condiciones de contorno: aplique la velocidad de entrada y velocidades cero dentro del ala.
El código
Así es como se ven las actualizaciones matemáticas principales en Python (vectorizadas para mejorar el rendimiento).
Paso A: Construir el término de la fuente de presión Esto representa el lado derecho (RHS) de la ecuación de Poisson según las velocidades actuales.
# b es el término fuente # u y v son matrices de velocidad actuales b[1:-1, 1:-1] = (rho * ( 1 / dt * ((u[1:-1, 2:] -tú[1:-1, 0:-2]) / (2 * dx) + (v[2:, 1:-1] -v[0:-2, 1:-1]) / (2 * dy)) – ((u[1:-1, 2:] -tú[1:-1, 0:-2]) / (2 * dx))**2 – 2 * ((u[2:, 1:-1] -tú[0:-2, 1:-1]) / (2 * dy) * (v[1:-1, 2:] -v[1:-1, 0:-2]) / (2 * dx)) – ((v[2:, 1:-1] -v[0:-2, 1:-1]) / (2 * dy))**2 ))
Paso B: Resolviendo la presión (iteración Jacobi) Repetimos para suavizar el campo de presión hasta que equilibre el término fuente.
para _ en rango (nit): pn = p.copy() p[1:-1, 1:-1] = ((pn[1:-1, 2:] +pn[1:-1, 0:-2]) * dy**2 + (pn[2:, 1:-1] +pn[0:-2, 1:-1]) * dx**2 – b[1:-1, 1:-1] * dx**2 * dy**2 ) / (2 * (dx**2 + dy**2)) # Condiciones de contorno: p=0 en los bordes (presión manométrica) p[:, -1] = 0; pag[:, 0] = 0; pag[-1, :] = 0; pag[0, :] = 0
Paso C: Actualización de la velocidad Finalmente, actualizamos la velocidad usando las ecuaciones de momento discretizadas explícitas.
un = u.copy() vn = v.copy() # Actualizar u (velocidad x) u[1:-1, 1:-1] = (un[1:-1, 1:-1] – ONU[1:-1, 1:-1] * dt / dx * (un[1:-1, 1:-1] – ONU[1:-1, 0:-2]) – vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] – ONU[0:-2, 1:-1]) – dt / (2 * rho * dx) * (p[1:-1, 2:] – pag[1:-1, 0:-2]) + nu * (dt / dx**2 * (un[1:-1, 2:] – 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) + dt / dy**2 * (un[2:, 1:-1] – 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]))) # Actualización v (velocidad y) v[1:-1, 1:-1] = (vn[1:-1, 1:-1] – ONU[1:-1, 1:-1] *dt/dx* (vn[1:-1, 1:-1] -vn[1:-1, 0:-2]) – vn[1:-1, 1:-1] *dt/dy* (vn[1:-1, 1:-1] -vn[0:-2, 1:-1]) – dt / (2 * rho * dy) * (p[2:, 1:-1] – pag[0:-2, 1:-1]) + nu * (dt / dx**2 * (vn[1:-1, 2:] – 2*vn[1:-1, 1:-1] +vn[1:-1, 0:-2]) + dt / dy**2 * (vn[2:, 1:-1] – 2*vn[1:-1, 1:-1] +vn[0:-2, 1:-1])))
Resultados: ¿vuela?
Ejecutamos este solucionador en un perfil de ala rígido con una entrada constante de campo lejano.
Observaciones cualitativas Los resultados se alinean con las expectativas físicas. Las simulaciones muestran alta presión debajo del ala y baja presión encima, que es exactamente el mecanismo que genera sustentación. Los vectores de velocidad muestran el flujo de aire acelerándose sobre la superficie superior (principio de Bernoulli).
Fuerzas: sustentación versus resistencia Al integrar el campo de presión sobre la superficie del ala, podemos calcular la sustentación.
El solucionador demuestra que las fuerzas de presión dominan las fuerzas de fricción viscosas por un factor de casi 1000 veces en el aire. A medida que aumenta el ángulo de ataque (de 0∘ a −20∘), la relación elevación-resistencia aumenta, coincidiendo con las tendencias observadas en túneles de viento y paquetes CFD profesionales como OpenFOAM.
Limitaciones y próximos pasos
Si bien crear este solucionador fue excelente para aprender, la herramienta en sí tiene sus limitaciones:
Resolución: las simulaciones 3D en una cuadrícula cartesiana son costosas desde el punto de vista computacional y requieren cuadrículas gruesas, lo que hace que los resultados cuantitativos sean menos confiables. Turbulencia: El solucionador es laminar; carece de un modelo de turbulencia (como k − ϵ ) necesario para flujos complejos o de alta velocidad. Difusión: Los esquemas de diferenciación a barlovento son estables pero numéricamente difusos, lo que potencialmente “borra” los detalles finos del flujo.
¿A dónde ir desde aquí? Este proyecto sirve como punto de partida. Las mejoras futuras podrían incluir la implementación de esquemas de advección de orden superior (como WENO), agregar modelado de turbulencia o pasar a métodos de volumen finito (como OpenFOAM) para un mejor manejo de malla alrededor de geometrías complejas. Existen muchas técnicas inteligentes para sortear la gran cantidad de escenarios que quizás desee implementar. ¡Éste es sólo el primer paso hacia una verdadera comprensión del CFD!