INSTITUTO POLITÉCNICO NACIONAL
“SIMULACIÓN DE NAVIER-STOKES EN
GEOMETRÍA VARIABLE”
-ROMERO MARTÍNEZ ADRIÁN YAHIR 7SV1
RESUMEN:
Este trabajo presenta la simulación numérica transitoria de un flujo viscoso e incompresible alrededor
de un obstáculo fijo dentro de un canal, utilizando el método de Diferencias Finitas implementado en
Python. El objetivo principal es estudiar la evolución temporal y espacial del campo de velocidad del
fluido (densidad ρ=1204.0 kg/m3, viscosidad cinemática ν=0.002 m2/s) para un Número de Reynolds
(Re=800). Este valor de Re sugiere un régimen de flujo laminar o de transición.
La simulación se lleva a cabo en un dominio de .6 m× .2 m con una malla discreta de 100×40 puntos.
Se emplea la solución de las Ecuaciones de Navier-Stokes en su forma no estacional, donde la velocidad
de entrada se calcula a partir del Número de Reynolds (u_entrada=Re⋅ν/L). El obstáculo (una placa
rectangular) se ubica en la esquina inferior izquierda. Se utilizan esquemas de diferencias finitas
centradas para las derivadas espaciales (términos convectivos, de presión y difusivos) y un esquema
explícito para la integración temporal con un paso de tiempo Δt=0.001 seg.
Los resultados se visualizan a través de mapas de contorno del campo de velocidad y perfiles de
velocidad en tres posiciones clave (antes, durante y después de la placa). El análisis se centra en la
distorsión del perfil de flujo y el desarrollo de la capa límite a lo largo del canal, cuantificando la
velocidad máxima, promedio y el espesor de la capa límite en cada punto de referencia.
MARCO TEORICO
El marco teórico de este proyecto se cimienta en la Mecánica de Fluidos y la Dinámica de Fluidos
Computacional (CFD).
1. Mecánica de Fluidos: Flujo Viscoso e Incompresible
El comportamiento físico del fluido está regido por dos principios fundamentales:
1.1. Ecuaciones de Navier-Stokes (ENS)
Estas ecuaciones representan la aplicación de las leyes de conservación de la masa y del momento
(Segunda Ley de Newton) a un volumen de control en el fluido. Para un flujo incompresible y
Newtoniano (donde la viscosidad es constante), las ENS se expresan de la siguiente forma vectorial:
• Conservación del Momento (Momentum):
∂t∂u+(u⋅∇)u=−ρ1∇p+ν∇2u+f
Donde:
• ∂t∂u: Término de acumulación (transitorio).
• (u⋅∇)u: Término de convección (transporte por el propio flujo).
• −ρ1∇p: Término de presión.
• ν∇2u: Término de difusión viscosa.
• u: Vector de velocidad. ν: Viscosidad cinemática. ρ: Densidad. p: Presión.
• Conservación de la Masa (Continuidad):
∇⋅u=0
Esta ecuación, al ser ∇⋅u=0, confirma la condición de incompresibilidad.
1.2. Número de Reynolds (Re)
El Número de Reynolds es un parámetro adimensional crucial que define el régimen del flujo:
Re=Fuerzas ViscosasFuerzas de Inercia=μρUL=νUL
Donde U es una velocidad característica y L es una longitud característica. Un Re=800 (utilizado en el
código) típicamente indica un flujo laminar estable o ligeramente inestable, donde las fuerzas viscosas
aún dominan sobre las turbulentas.
1.3. Capa Límite y Separación del Flujo
La interacción del fluido con la superficie sólida del obstáculo genera una capa límite, una delgada
región adyacente a la pared donde los efectos de la viscosidad son significativos y la velocidad del fluido
pasa de cero (condición de no deslizamiento) a la velocidad de la corriente libre. La presencia del
obstáculo y los gradientes de presión pueden provocar el fenómeno de separación del flujo (la capa
límite se separa de la pared), creando remolinos o estelas que son cruciales para el análisis de los
perfiles de velocidad.
2. Dinámica de Fluidos Computacional (CFD)
El código implementa una solución numérica de las ENS:
2.1. Método de Diferencias Finitas (FDM)
El FDM es el corazón de la simulación. Consiste en:
• Discretización del Dominio: Dividir el dominio continuo (L×H) en una rejilla o malla discreta
(Δx, Δy).
• Discretización de las Ecuaciones: Reemplazar las derivadas parciales de las ENS con
aproximaciones de diferencias finitas (diferencias centradas de segundo orden, como se observa
en el código). Por ejemplo, la derivada ∂x∂u se aproxima como 2Δxuj+1−uj−1.
2.2. Esquemas de Integración Temporal
El código utiliza una solución transitoria (dependiente del tiempo), implementando un esquema
explícito (posiblemente Forward Euler o similar) donde los valores futuros de u y v se calculan
directamente a partir de los valores en el paso de tiempo actual, como se ve en la estructura:
Unuevo=Uviejo+Δt⋅(Teˊrminos de la Ecuacioˊ n)
2.3. Condiciones de Contorno
Las condiciones de contorno son esenciales para definir el problema físico:
• Entrada : Velocidad horizontal uniforme (U=u_entrada) y velocidad vertical nula (V=0).
• Paredes (Superior/Inferior): Condición de no deslizamiento (U=0,V=0).
• Obstáculo: Condición de no deslizamiento (U=0,V=0).
• Salida : Se aplica una condición de salida de flujo completamente desarrollado o "cero
gradiente" (U_final=U_penultimo y V_final=V_penultimo), lo que permite que el flujo "salga"
del dominio sin reflejar perturbaciones.
• Condición de Presión: En este método, la presión se calcula en un paso separado (aunque el
cálculo no se muestra explícitamente en el solver de momento, los términos de presión se usan
para actualizar la velocidad). El método de proyección (que resuelve la ecuación de Poisson para
la presión) es el estándar para acoplar la velocidad y la presión en flujos incompresibles.
DISCRETIZACIÓN
U=(convective_u)=(u∂x∂u+v∂y∂u)=Ui,jo2ΔxUi,j+1o
−Ui,j−1o+Vi,jo2ΔyUi+1,jo−Ui−1,jo
V=(convective_v)=(u∂x∂v+v∂y∂v)=Ui,jo2ΔxVi,j+1o
−Vi,j−1o+Vi,jo2ΔyVi+1,jo−Vi−1,jo
DESARROLLO
Lo primero que hacemos es definir nuestro mundo. Usamos librerías como NumPy para todas las
matemáticas, porque el corazón de esto son las matrices, y Matplotlib para dibujar los resultados.
Definimos el tamaño de nuestra piscina (Largo L=6, Alto H=2) y las propiedades del fluido: su
viscosidad (ν) y su densidad (ρ). Luego, usamos el Número de Reynolds (Re=100) para calcular la
velocidad de entrada (u_entrada) que necesitamos, asegurándonos de que nuestro flujo se comporte
como queremos. También establecemos la cuadrícula, o malla, con 100 puntos a lo largo y 40 a lo alto
(nx, ny), lo que determina cuán detallada será la simulación. El paso de tiempo (Δt=0.01) es clave: es el
interruptor que nos permite avanzar en el tiempo sin que la simulación se vuelva caótica. Finalmente,
definimos el obstáculo como un bloque en la esquina inferior izquierda de la malla.
Las líneas nx, ny = 200, 200 y dx, dy = L/nx, H/ny establecen la malla o cuadrícula de cálculo, que es la
base de toda la simulación. Al usar 200 puntos tanto en la dirección horizontal (nx) como en la vertical
(ny), estamos dividiendo el dominio físico que definimos antes (L y H) en una gran cantidad de
pequeños cuadros, y el dx y dy calculan el tamaño exacto de esos cuadritos. Cuanto mayor sea este
número, más detallada será nuestra simulació[Link] líneas nx_en, ny_en = int(0.2 * nx), int(0.5 * ny)
definen el tamaño del obstáculo que el fluido debe rodear. Estamos diciendo que el obstáculo debe
ocupar el 20% de los puntos en x y el 50% de los puntos en y de nuestra malla, lo que determina su
ubicación y tamaño en las coordenadas de la cuadrícula, que por lo general está en la esquina de entrada
en este tipo de simulaciones.
Re = DESEADO Y PROPUESTO 100,300,800
L, H = 6, 2
v = 0.002
rho = 1204.0
u_entrada = Re * v / L
print(f"Velocidad de entrada calculada: u_ENTRADA = {u_entrada:.4f} m/s")
tiempo_total = 10.0
dt = 0.01
pasos_totales = int(tiempo_total / dt)
intervalo = 100
max_iter_poisson = 50
nx, ny = 200, 200
dx, dy = L/nx, H/ny
nx_en, ny_en = int(0.2 * nx), int(0.5 * ny)
Una vez que tenemos el mapa, tenemos que llenar cada punto de la cuadrícula con valores iniciales para
las tres variables principales que vamos a resolver: la velocidad y la presión. El código usa numpy para
crear estas grandes matrices. La línea U = [Link]((ny, nx), u_entrada) crea la matriz de
velocidad horizontal (U), y la inicializa completamente con la u_entrada que calculamos al inicio, como
si todo el fluido ya estuviera moviéndose a esa velocidad. Luego, V = [Link]((ny, nx)) crea la
matriz de velocidad vertical (V) y la inicializa a cero, asumiendo que el fluido inicialmente solo se
mueve horizontalmente. Finalmente, P = [Link]((ny, nx)) crea la matriz de presión (P) e
igual la inicializa a cero en todas partes. Estas matrices son cruciales porque serán rellenadas con los
resultados reales en cada paso de tiempo del simulador. Por último, las listas vacías
tiempos_guardados = [] y campos_guardados = [] son simplemente el lugar donde
guardaremos los resultados de la simulación.
U = [Link]((ny, nx), u_entrada)
V = [Link]((ny, nx))
P = [Link]((ny, nx))
tiempos_guardados = []
campos_guardados = []
Esta parte del código es absolutamente vital porque aquí definimos las reglas de interacción en los
límites de nuestra simulación; son las condiciones de frontera que le dicen al fluido y a la presión cómo
deben comportarse cuando tocan una pared o un límite. En la función set_CONDICIONES, que maneja
las velocidades (U y V), decimos que en las paredes superior e inferior, y en todas las celdas del
obstáculo, la velocidad es cero (U[0, :] = 0; V[0, :] = 0), lo que representa la condición de no
deslizamiento, que significa que el fluido se pega a las superficies sólidas. En la entrada (borde
izquierdo, índice 0), la velocidad horizontal se fija a u_entrada y la vertical a cero, asegurando un flujo
constante que entra al dominio. En la salida (borde derecho, índice -1), le decimos que la velocidad
debe ser igual a la de la fila adyacente (U[:, -1] = U[:, -2]), lo que permite al fluido salir libremente sin
ser perturbado. Por otro lado, la función set_CONDICIONES_PRESION maneja el campo de presión
(P), y para todos los bordes, incluyendo el obstáculo, imponemos una condición de gradiente cero (P[:,
0] = P[:, 1]), lo que permite que la presión en el borde se ajuste al valor de la celda vecina.
def set_CONDICIONES(U, V, u_entrada):
U[0, :] = 0; U[-1, :] = 0
V[0, :] = 0; V[-1, :] = 0
U[:, 0] = u_entrada
V[:, 0] = 0
U[:, -1] = U[:, -2]
V[:, -1] = V[:, -2]
U[:ny_en, :nx_en] = 0
V[:ny_en, :nx_en] = 0
return U, V
def set_CONDICIONES_PRESION(P):
P[:, 0] = P[:, 1]
P[:, -1] = P[:, -2]
P[0, :] = P[1, :]
P[-1, :] = P[-2, :]
P[:ny_en, :nx_en] = P[ny_en, nx_en]
return P
Esta función resuelve la Ecuación de Poisson de forma iterativa para encontrar el campo de presión que
garantiza la incompresibilidad. Primero, calcula la divergencia de las velocidades intermedias (D)
usando diferencias centrales, que se convierte en el término fuente. Luego, en el ciclo de iteración,
aproxima el Laplaciano de la presión (P_xx + P_yy) a partir de los valores de los vecinos.
Finalmente, utiliza estos términos y la divergencia D para actualizar la presión (Pn), forzando que la
divergencia se anule, y repite el proceso hasta alcanzar la convergencia o el número máximo de
iteraciones.
La función SOLVER es el motor principal que avanza la simulación paso a paso. En cada step, primero
calcula una velocidad intermedia (Ue) considerando la convección (inercia) y la difusión (viscosidad),
sin incluir la presión. Luego, utiliza la función POISSON para encontrar el campo de presión necesario
para que el fluido sea incompresible. Finalmente, corrige las velocidades intermedias Ue y Ve usando el
gradiente de presión (∇P) recién calculado para obtener las velocidades finales U y V para el siguiente
paso, cerrando el ciclo.
def POISSON(P, U_e, V_e, rho, dt, dx, dy, max_iter):
ny, nx = U_e.shape
Pn = [Link]()
for iter_count in range(max_iter):
Po = [Link]()
D = [Link]((ny-2, nx-2))
for i in range(1, ny-1):
for j in range(1, nx-1):
dU_dx = (U_e[i, j+1] - U_e[i, j-1]) / (2 * dx)
dV_dy = (V_e[i+1, j] - V_e[i-1, j]) / (2 * dy)
D[i-1, j-1] = dU_dx + dV_dy
for i in range(1, ny - 1):
for j in range(1, nx - 1):
if i < ny_en and j < nx_en:
continue
P_xx = (Po[i, j+1] + Po[i, j-1]) / (dx**2)
P_yy = (Po[i+1, j] + Po[i-1, j]) / (dy**2)
Pn[i, j] = (P_xx + P_yy - D[i-1, j-1] * rho / dt) / (2 /
dx**2 + 2 / dy**2)
Pn = set_CONDICIONES_PRESION(Pn)
return Pn
Una vez que el ciclo de SOLVER termina, toda esta sección se encarga de guardar los resultados y
mostrarlos. Primero, el bloque if step % save_interval == 0: se ejecuta periódicamente
para guardar el estado de los campos de velocidad (U,V) y presión (P) en las listas
campos_guardados y tiempos_guardados, y nos da una actualización impresa del
[Link]és de que la función SOLVER ha corrido durante todo el tiempo_total, se
recuperan los resultados finales en U_final, V_final, y P_final. La última sección prepara el
gráfico: primero crea los ejes de coordenadas X e Y con [Link], y luego usa Matplotlib para
visualizar el flujo. El comando [Link] dibuja el mapa de colores de la velocidad final, y
[Link] agrega las flechas de vectores para mostrar la dirección y magnitud del flujo.
if step % save_interval == 0:
tiempo_actual = step * dt
tiempos_guardados.append(tiempo_actual)
campos_guardados.append({'U': [Link](), 'V': [Link](), 'P':
[Link](), 'tiempo': tiempo_actual})
print(f"Paso {step:4d}, Tiempo: {tiempo_actual:6.2f} s, U_max:
{[Link](U):.4f} m/s, P_max: {[Link](P):.2f} Pa")
return U, V, P, tiempos_guardados, campos_guardados
U_final, V_final, P_final, tiempos, campos = SOLVER(U, V, P, u_entrada,dt,
pasos_totales, intervalo)
x = [Link](0, L, nx)
y = [Link](0, H, ny)
X, Y = [Link](x, y)
[Link](figsize=(10, 8))
contour = [Link](X, Y, U_final, levels=50, cmap='viridis')
[Link](X,Y,U_final,V_final)
[Link](contour, label='Presión [Pa]')
[Link](f'Campo de Velocidad (t = {tiempo_total} s)\nRe = {Re}')
[Link]('X [m]')
[Link]('Y [m]')
obstaculo_x = [0, x[nx_en], x[nx_en], 0, 0]
obstaculo_y = [0, 0, y[ny_en], y[ny_en], 0]
[Link](obstaculo_x, obstaculo_y, 'gray', alpha=0.7, label='Obstáculo')
[Link]()
import numpy as np
import [Link] as plt
Re = 800
L, H = 6, 2
v = 0.002
rho = 1204.0
u_entrada = Re * v / L
print(f"Velocidad de entrada calculada: u_ENTRADA = {u_entrada:.4f} m/s")
tiempo_total = 10.0
dt = 0.01
pasos_totales = int(tiempo_total / dt)
intervalo = 100
max_iter_poisson = 50
nx, ny = 200, 200
dx, dy = L/nx, H/ny
nx_en, ny_en = int(0.2 * nx), int(0.5 * ny)
U = [Link]((ny, nx), u_entrada)
V = [Link]((ny, nx))
P = [Link]((ny, nx))
tiempos_guardados = []
campos_guardados = []
def set_CONDICIONES(U, V, u_entrada):
U[0, :] = 0; U[-1, :] = 0
V[0, :] = 0; V[-1, :] = 0
U[:, 0] = u_entrada
V[:, 0] = 0
U[:, -1] = U[:, -2]
V[:, -1] = V[:, -2]
U[:ny_en, :nx_en] = 0
V[:ny_en, :nx_en] = 0
return U, V
def set_CONDICIONES_PRESION(P):
P[:, 0] = P[:, 1]
P[:, -1] = P[:, -2]
P[0, :] = P[1, :]
P[-1, :] = P[-2, :]
P[:ny_en, :nx_en] = P[ny_en, nx_en]
return P
def POISSON(P, U_e, V_e, rho, dt, dx, dy, max_iter):
ny, nx = U_e.shape
Pn = [Link]()
for iter_count in range(max_iter):
Po = [Link]()
D = [Link]((ny-2, nx-2))
for i in range(1, ny-1):
for j in range(1, nx-1):
dU_dx = (U_e[i, j+1] - U_e[i, j-1]) / (2 * dx)
dV_dy = (V_e[i+1, j] - V_e[i-1, j]) / (2 * dy)
D[i-1, j-1] = dU_dx + dV_dy
for i in range(1, ny - 1):
for j in range(1, nx - 1):
if i < ny_en and j < nx_en:
continue
P_xx = (Po[i, j+1] + Po[i, j-1]) / (dx**2)
P_yy = (Po[i+1, j] + Po[i-1, j]) / (dy**2)
Pn[i, j] = (P_xx + P_yy - D[i-1, j-1] * rho / dt) / (2 /
dx**2 + 2 / dy**2)
Pn = set_CONDICIONES_PRESION(Pn)
return Pn
def SOLVER(U, V, P, u_entrada, dt, total_steps, save_interval):
ny, nx = [Link]
for step in range(total_steps):
Uo, Vo = [Link](), [Link]()
Ue = [Link]()
Ve = [Link]()
for i in range(1, ny - 1):
for j in range(1, nx - 1):
if i < ny_en and j < nx_en: continue
convective_u = Uo[i, j] * (Uo[i, j+1] - Uo[i, j-1]) /
(2*dx) + \
Vo[i, j] * (Uo[i+1, j] - Uo[i-1, j]) /
(2*dy)
convective_v = Uo[i, j] * (Vo[i, j+1] - Vo[i, j-1]) /
(2*dx) + \
Vo[i, j] * (Vo[i+1, j] - Vo[i-1, j]) /
(2*dy)
diffusive_u = v * ( (Uo[i, j+1] - 2*Uo[i, j] + Uo[i, j-1])
/ (dx**2) +
(Uo[i+1, j] - 2*Uo[i, j] + Uo[i-1, j])
/ (dy**2) )
diffusive_v = v * ( (Vo[i, j+1] - 2*Vo[i, j] + Vo[i, j-1])
/ (dx**2) +
(Vo[i+1, j] - 2*Vo[i, j] + Vo[i-1, j])
/ (dy**2) )
Ue[i, j] = Uo[i, j] + dt * (-convective_u + diffusive_u)
Ve[i, j] = Vo[i, j] + dt * (-convective_v + diffusive_v)
Ue, Ve = set_CONDICIONES(Ue, Ve, u_entrada)
P = POISSON(P, Ue, Ve, rho, dt, dx, dy, max_iter_poisson)
for i in range(1, ny - 1):
for j in range(1, nx - 1):
if i < ny_en and j < nx_en: continue
grad_P_x = (P[i, j+1] - P[i, j-1]) / (2*dx)
U[i, j] = Ue[i, j] - (dt / rho) * grad_P_x
grad_P_y = (P[i+1, j] - P[i-1, j]) / (2*dy)
V[i, j] = Ve[i, j] - (dt / rho) * grad_P_y
U, V = set_CONDICIONES(U, V, u_entrada)
if step % save_interval == 0:
tiempo_actual = step * dt
tiempos_guardados.append(tiempo_actual)
campos_guardados.append({'U': [Link](), 'V': [Link](), 'P':
[Link](), 'tiempo': tiempo_actual})
print(f"Paso {step:4d}, Tiempo: {tiempo_actual:6.2f} s, U_max:
{[Link](U):.4f} m/s, P_max: {[Link](P):.2f} Pa")
return U, V, P, tiempos_guardados, campos_guardados
U_final, V_final, P_final, tiempos, campos = SOLVER(
U, V, P, u_entrada, dt, pasos_totales, intervalo
)
x = [Link](0, L, nx)
y = [Link](0, H, ny)
X, Y = [Link](x, y)
[Link](figsize=(10, 8))
contour = [Link](X, Y, U_final, levels=50, cmap='viridis')
[Link](X,Y,U_final,V_final)
[Link](contour, label='Presión [Pa]')
[Link](f'Campo de Velocidad (t = {tiempo_total} s)\nRe = {Re}')
[Link]('X [m]')
[Link]('Y [m]')
obstaculo_x = [0, x[nx_en], x[nx_en], 0, 0]
obstaculo_y = [0, 0, y[ny_en], y[ny_en], 0]
[Link](obstaculo_x, obstaculo_y, 'gray', alpha=0.7, label='Obstáculo')
[Link]()
GRAFICOS: