2.3. Hidrodinámica estacionaria en dos dimensiones.

• Teoría: Medios contínuos, hidrodinámica.
• Ténica numérica: Resolución de ecuaciones en derivadas parciales de tipo elíptico por el método de la relajación de Gauss-Seidel.
• Programación: Fortran, C, DynamicLattice.

Introducción teórica:

Queremos conocer el comportamiento estacionario de un fluido incompresible (por ejemplo el agua) cuando se encuentra en su camino con un obstáculo. Nos interesa saber el tipo de regímenes que tiene el fluido dependiendo de su flujo incidente y la presión que se ejerce sobre el objeto con el fin de diseñar geometrías que minimicen dicha presión.

Supondremos que el objeto es invariante por traslaciones en una de las direcciones espaciales (como por ejemplo z). De esta forma el movimiento del fluido en el plano (x,y) no será trivial.

El estado de un fluido queda determinado al conocer en cada punto (x,y) y para cualquier instante t su densidad, r(x,y;t), su vector velocidad v(x,y;t) y su presión P(x,y;t). Este conjunto de cuatro variables cumplen las siguientes ecuaciones diferenciales en derivadas parciales:

donde D simboliza el operador gradiente sobre las coordenadas espaciales. La primera de las ecuaciones se denomina ecuación de continuidad y refleja que la variación en el tiempo de la masa en un punto del espacio solo puede cambiar debido a una corriente local de masa (conservando globalmente la misma).

La segunda ecuación es la famosa ecuación de Navier-Stokes que refleja la conservación del momento lineal. Esto es, la velocidad puede cambiar en el tiempo debido a (i) convección (primer término), (ii) a variaciones locales de la presión (segundo término) y (iii) a la acción de fuerzas viscosas (último término). n es la denominada viscosidad cinemática y se supone constante.

Las anteriores dos ecuaciones son insuficientes para obtener una solución pues nos falta una ecuación que relacione la presión y la densidad (ecuación de estado). Normalmente las ecuaciones de estado dependen de la temperatura, por lo que si permitimos variaciones locales de la temperatura deberíamos incluir también una ecuacion de continuidad para la energía.

En nuestro modelo supondremos: (1) n = cte; (2) la temperatura es constante (no hay flujos de energía); (3) la densidad es constante (fluido incompresible) y (4) estado estacionario i.e. d_tr = 0 y d_tv = 0. De esta forma, las ecuaciones en derivadas parciales que hemos de resolver son:

Si definimos v(x,y) = (u(x,y),v(x,y)) podemos reescribir las anteriores ecuaciones como:

donde las incógnitas son u, v y P. Podemos simplificar bastante las ecuaciones si introducimos dos nuevas funciones:

• (1) Función corriente (f): Es aquella función tal que

  u = dy f    v = -dx f

  (8)

  Notar que esta función f existe y es única salvo una constante. La ventaja de su utilización es que la primera de las ecuaciones se cumple idénticamente cuando la expresamos en términos de f. f tiene sentido físico: el vector velocidad v(x,y) es tangente a la superficie f(x,y) = cte. Esta propiedad es sencilla de demostrar: es obvio que se cumple v·Df = 0, además, si nos movemos un dr sobre la superficie f(x,y) = cte, el valor de f no variará, esto es, Df·dr = 0. Asi pues, el gradiente de f es perpendicular a la superficie f = cte y por otra parte es perpendicular al vector velocidad, por lo que éste último debe de ser tangente a la superficie f(x,y) = cte.

• (2) Función vorticidad (y): es aquella que nos da idea de la existencia de vórtices en el fluido:

  y = -D×v = dy u-dx v

  (9)

Utilizando estas dos definiciones podemos reescribir las ecuaciones de Navier-Stokes de la siguiente forma:

La demostración de la primera equation no es más que substituir la relación entre las derivadas de f y las velocidades. La segunda se obtiene tomando las dos ecuaciones de Navier-Stokes y derivando la primera respecto d_y y la segunda respecto d_x y restando las ecuaciones resultantes entre si. Finalmente, la tercera ecuación se obtiene derivando la primera respecto d_x, la segunda respecto d_y y sumando las ecuaciones. En todos los casos se utiliza la propiedad d_x u = d_y v. Notar que las tres ecuaciones resultantes son tales que podemos resolver únicamente las dos primeras para conocer el campo de velocidades. Si queremos conocer la presión, tendríamos que resolver posteriormente la tercera.

Problemas:

• Obligatorio: Obtener los campos f y y solución de las ecuaciones hidrodinámicas para un fluido incompresible correpondiente al caso de un flujo sorteando un obstáculo rectangular perpendicular a su velocidad incidente.
• Voluntario: Estudiar el efecto de la velocidad del flujo incidente en el comportamiento de los campos f y y. Obtener la velocidad crítica que separe los regímenes laminar y turbulento.

Método numérico:

El espacio correspondiente al movimiento del fluido en presencia del obstáculo lo definimos por un retículo bidimensional cuyos nodos tienen coordenadas: x_i = ih donde i = 1,..,N_x e y_j = jh donde j = 1,..,N_y. Donde h es el parámetro que define la discretización del espacio. Esta discretización la reflejamos en las derivadas parciales utilizando la siguiente identificación:

donde f_i,j = f(ih,jh). Nos vamos a concentrar en la resolución de los campos corriente y vorticidad. Utilizando las anteriores dicretizaciones, las ecuaciones que hemos de resolver son:

Es conveniente en este momento introducir un conjunto de nuevas variables reescaladas. Para elegirlas nos fijamos en la relación dimensional entre f y u: [u] = [f]/[y]. De esta forma, si lo más natural es medir las distancias en unidades h: y = y^(s) h y las velocidades en unidades del flujo incidente: u = u^(s)V₀, obtenemos que [u^(s)] = [f]/[y^(s)][V₀h] asi, la forma de obtener una f escalada que sea adimensional es definir f = f^(s) V₀ h. De igual forma razonamos el comportamiento de f, sabemos que [y] = [u]/[y], asi, sustituyendo las velocidades y posiciones reescaladas obtenemos [y] = [u^(s)][V₀]/[y^(s)][h] de forma que y = y^(s) V₀/h. Se puede comprobar (por completitud) que P = P^(s)rV₀². En función de las nuevas variables, las ecuaciones discretas quedan:

donde el único parámetro que queda es R = V₀ h/n. Notar que antes de realizar la discretización uno podía haber definido un parámetro físico adimensional, el llamada número de Reynolds: Re = 2WV₀/n. Donde 2W es la longitud transversal del obstáculo con respecto al flujo. Notar que Re = 2WR/h, esto es, la fenomenología macroscópica está caracterizada por Re y W, de forma que si variamos h hemos de variar R para describir lo mismo. Asi los resultados correctos los obtendremos formalmente como la extrapolación de las tendencias de diferentes simulaciones con distintos reticulos y sus correspondientes R.

En general, la forma de resolver las anteriores ecuaciones discretas es mediante un proceso de relajación iterativa. Ilustrémoslo con un ejemplo. Supongamos que hemos de resolver la ecuación discretizada:

Si despejamos f_i obtenemos:

esto nos da la idea para construir una nueva configuración de f en función de la antigua:

donde se introduce el parámetro w para moderar la evolución de la iteración sobre todo en las fronteras donde inicialmente las variaciones entre f_i y las de sus vecinos próximos son demasiado acusadas.

La aplicación de esta técnica relajacional en nuestras dos ecuaciones acopladas es sencillo. Dadas las configuraciones de f^(s) y y^(s), primero realizaremos una iteración relajacional sobre las f^(s) suponiendo y^(s) fijas y posteriormente haremos una iteración sobre y^(s) manteniendo las f^(s) fijas. Iremos iterando estos dos pasos hasta que logremos una configuración estable y estacionaria bajo las dos iteraciones.

Antes de hacer explícito el algoritmo de relajación, hemos de discutir con cierta profundidad las condiciones de contorno que hemos de aplicar a nuestro problema (Ver figura). Hemos de tener en cuenta varias ideas:

• (1) Las fronteras que tenemos son A,B,...,H, donde el eje X es un eje de simetría bajo reflexión de nuestro problema.
• (2) El fluido se supone que incide a través de F y se escapa por H.
• (3) G se supone lo suficientemente alejado del obstáculo para que el fluido a partir de alli se comporte como si no existiese el obstáculo.
• (4) La velocidad normal a las superficies del obstáculo es cero.
• (5) Las condiciones de contorno sobre el obstaculo son las denominadas no slip boundary conditions que consisten en que las velocidades normales y transversales son cero en el obstáculo.

Con las anteriores ideas en mente, podemos construir las condiciones de contorno para la función corriente y vorticidad. Comencemos con la primera.

Condiciones de contorno para la función corriente (f):

• En A: Puesto que A es un eje de simetría de nuestro problema, la componente y de la velocidad en cada punto de A debe de ser cero. Asi d_xf = 0, lo que implica f(x,y) = f(y). Si tomamos como origen de coordenadas la intersección entre E y F f = f(0) = cte. Esa constante es arbitraria y señala el cero del campo corriente, el fijarlo a cero no influye en el comportamiento de las velociades puesto que éstas son derivadas de f. En conclusión f = 0 y por lo tanto f^(s)_i,j = 0 en A.
• En B: La velocidad normal al obstáculo ha de ser cero, luego u = 0 lo que implica d_yf = 0. Asi f(x,y) = g(x), en particular como x es constante en B, f será constante en B. Su valor coincidirá con el valor de f en A puesto que tienen un punto en común. Asi f = 0 y por lo tanto f^(s)_i,j = 0 en B.
• En C: Se argumenta lo mismo que en B solo que ahora v = 0. Asi f = 0 y por lo tanto f^(s)_i,j = 0 en C.
• En D: Lo mismo que en B. Asi f = 0 y por lo tanto f^(s)_i,j = 0 en B.
• En E: Lo mismo que en A. Asi f = 0 y por lo tanto f^(s)_i,j = 0 en A.
• En F: Se supone que en F el fluido se comporta como si estuviesemos a distancia infinita del obstáculo, esto es, el fluido se mueve horizontalmente con velocidad constante. Asi d f/d x = 0 en F. Lo que implica que f^(s)_0,j = f^(s)_1,j. Notar que esta condición de contorno se diferencia de las anteriores. La estrategia a seguir será la de realizar una iteración sobre f^(s) dadas unas condiciones de contorno y al acabar la iteración cambiar la condición de contorno en F. El valor inicial que se le puede dar a f^(s) es el de flujo libre: f^(s)_0,j = j.
• En G: Al igual que en F, se supone que el fluido tiene velocidad u = V₀ y v = 0 sobre G. De esta forma d_yf = V₀ lo que implica que f = V₀y. Asi pues: f^(s)_i,Ny+1 = f^(s)_i,Ny en G y utilizaremos la misma estrategia que en F en la actualización de su valor, tomando como condición inicial f^(s)_i,Ny+1 = N_y+1.
• En H: Igual que en F solo que ahora f^(s)_Nx+1,j = f^(s)_Nx,j.

Notar que no necesitamos utilizar el hecho de que la velocidad tangencial sobre el obstáculo es cero. Esta información la utilizares en las condiciones de contorno para y.

Condiciones de contorno para la función corriente (y):

• En A: Por simetría, v en A es cero, por lo que d_x v = 0 en A. También, en las cercanías de y = 0, u(x,y) debe de ser par en y lo que implica d_y u = 0 en A. Asi y = 0 en A y por lo tanto y^(s) = 0 en A.

• En B: Sabemos que en las cercanías del obstáculo y para una j dada se cumple:

  f((i+1)h,jh) = f(ih,jh)+h df d x |ih,jh+ h2 2 d2f d x2 + ...

  (22)

  pero como v = 0 sobre el obstáculo tenemos que d_xf|_ih,jh = 0. Por otra parte, sabemos que d_xxf = -d_x v = y puesto que u es constante (e igual a cero) sobre B. Asi pues, despejando la segunda derivada en función de los incrementos de f, obtenemos el valor de y sobre el obstáculo B:

  y(s)i,j = 2f(s)i+1,j

  (23)

  Notar otra vez que para definir las condiciones de contorno sobre B, hemos de conocer primero f. Donde hemos utilizado el hecho de que f sobre el obstáculo es cero.

• En C: Se razona de forma similar a B y obtenemos que:

  y(s)i,j = 2f(s)i,j+1

  (24)

• En D: Igual que en B pero sustituyendo i+1 por i-1.

  y(s)i,j = 2f(s)i-1,j

  (25)

• En E: Igual que en A, luego y^(s) = 0 en E.
• En F: Sabemos que u no depende de y y v es igual a cero (lo mismo que su primera derivada) asi, y^(s) = 0 en F.
• En G: Igual que en F.
• En H: Tomamos como condición de frontera el que y no varíe en x, i.e. d_x y = 0, que no es más que la misma condición de contorno que f en esa pared. Asi: y^(s)_Nx+1,j = y^(s)_Nx,j. Al igual que en el caso de f, actualizaremos esta condición de contorno una vez calculada una iteración en y.

Asi pues, el algoritmo para resolver este problema es:

• (0) Dar valores iniciales de fluido libre.

  (a) A f^(s) en la red: f^(s)_i,j = j y sus condiciones de contorno: f^(s)_A = f^(s)_B = f^(s)_C = f^(s)_D = f^(s)_E = 0, f^(s)_i,j = j en F, G y H.

  (b) A y^(s) en la red y en sus condiciones de contorno : y^(s)_i,j = 0.

• (1) Fijar los valores w₁ y w₂ para las iteraciones de relajación de f y y respectivamente, asi como R.
• (2) Realizar una iteración de relajación de f^(s).

  (a) Contruir una nueva matriz

  fi,j = (1-w1)f(s)i,j+ w1 4 [f(s)i+1,j +f(s)i-1,j +f(s)i,j+1 +f(s)i,j-1-y(s)i,j ]

  (26)

  donde (i,j) pertencen al interior de la red (no son puntos frontera).

  (b) Substituir los antiguos valores de f por los nuevos: f^(s)_i,j = f_i,j.

  (c) Actualizar las condiciones de contorno F,G y H de f^(s): f^(s)_0,j = f^(s)_1,j, f^(s)_i,Ny+1 = f^(s)_i,Ny+1 y f^(s)_Nx+1,j = f^(s)_Nx,j.

• (3) Realizar una iteración de relajación de y^(s).

  (a) Modificar las condiciones de contorno en B: y^(s)_i,j = 2f^(s)_i+1,j, en C: y^(s)_i,j = 2f^(s)_i,j+1 y en D: y^(s)_i,j = 2f^(s)_i-1,j.

  (b) Construir una nueva matriz

  gi,j = (1-w2)y(s)i,j+ w2 4 [y(s)i+1,j+y(s)i-1,j+y(s)i,j+1+y(s)i,j-1

  (27)

  - R 4 [ (f(s)i,j+1-f(s)i,j-1)(y(s)i+1,j-f(s)i-1,j)

  (28)

  -(f(s)i+1,j-f(s)i-1,j)(y(s)i,j+1-f(s)i,j-1)]]

  (29)

  (c) Substituir los antiguos valores de y por los nuevos: y^(s)_i,j = g_i,j.

  (d) Actualizar la condición de contorno en H de y^(s): y^(s)_Nx+1,j = y^(s)_Nx,j.

• (4) Ir a (2).

Ayudas:

• Comprobar que si disminuimos el tamaño del retículo pero conservamos los valores físicos del problema, obtenemos lo mismo.
• Para simular números de Reynolds grandes, ir incrementado progresivamente los números de Reynolds y tomar como condición inicial la configuración estacionaria conrrespondiente a un valor algo más pequeño.